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1.  INTRODUCTION 


In  preparation  for  analysis  of  the  CRRES  databases  (to  be  taken  from  1990  to  1995)  and,  in 
particular,  for  an  updated  static  radiation  belt  model,  a  comparison  has  been  made  of  four 
magnctospheric  internal  magnetic  field  models.  Since  there  is  interest  in  radiation  belt 
models  out  to  geosynchronous  altitudes,  it  will  be  necessary  to  consider  the  external  magnetic 
field  models  as  well,  but  for  now  just  the  internal  models  will  be  discussed.  The  four  models 
reviewed  here  are  Barradough,  MAGSAT,  1GRF  1965,  and  Cain.  Each  of  these  has  an 
optimal  time  period  associated  with  it:  Barradough  1975,  MAGSAT  1980,  IGRF  1985,  and 
Cain  1990. 


In  the  subsequent  sections  of  this  report,  the  derivations  of  these  models  will  be  summarized. 
They  are  all  based  on  the  spherical  harmonic  expansion  of  the  magnetic  field’s  scalar 
potential.  The  data  used  to  determine  each  set  of  confidents  will  be  briefly  discussed.  Then 
the  models  will  be  evaluated  by  comparing  the  magnetic  field  and  L-shell  values  obtained 
from  each  of  the  models.  Finally,  various  computer  codes  which  use  these  models  to 
evaluate  the  field  and  L-shell,  will  be  compared  and  evaluated  in  terms  of  their  performance 
(i.e.  accuracy  and  speed). 


'~A  primary  concern  for  radiation  belt  models  is  reduction  of  the  number  of  binning 
parameters  by  using  B-L  coordinates, [McILwain,  1961].  Thus,  these  four  models  were 
compared  on  the  basis  of  the  total  magnetic  field,  its  components,  and  L-shell  values 
determined  from  each  model’s  unique  set  of  spherical  harmonic  coefficients.  The 
comparisons  were  made  to  high  latitudes  (±75°  in  geographic  coordinates),  over  all 
longitudes,  and  over  &  large  altitude  range  (0  km  -  40,000  km).  This  was  done  in  the  interest 
of  completeness.  However,  the  CRRESmrbits  will  fall  within  a  latitude  range  between  ±20°. 
Also,  beyond  approximately  15,000  km  the  magnetic  field  due  to  external  sources  becomes 
significant.  Thus,  the  focus  of  the  analysis  of  these  models  will  be  within  more  limited 
latitude  and  altitude  ranges.  _  .  *  ✓'  j  ,  r^l  '  ;  jf . 

V;-  ' 

In  addition  to  the  analysis  of  the  models  themselves,  ihe  software  in  which  they  are  used  h 
been  reviewed  as  well.  There  are  essentially  four  cooes  which  have  been  considered,  two  of  - 
which  are  currently  in  use  at  the  Geophysics  laboratory  (GL)  and  two  of  which  were  sent  6jr  Y> 
Dr.  J.  C  Cain.  The  two  codes  at  GL  are  OPTRACE  and  HMIN.  OPTRACE  is  a  more 
extensive  package  than  HMIN.  Initially,  OPTRACE  was  used  in  the  model  comparisons. 
However,  when  the  codes  themselves  were  to  be  evaluated,  a  modified,  version  of  HMIN  was 
used  to  give  a  more  accurate  comparison  of  the  CPU  requirements  i/f  each  of  the  codes. 
These  codes  trace  along  the  field  line  in  cartesian  corrdinates.  The  two  codes  sent  by  Cain 
were  reconstructed  from  work  done  by  Kluge  in  the  early  1970s.  One  of  these  integrates 
over  the  field  line  to  find  L  using  inverse  coordinates  centered  on  the  dipole.  This  package 
will  be  referred  to  as  Kluge/Cain  Integration.  The  other  package  is  a  short  routine  which 
uses  interpolation  tables  based  on  the  Cain  1990  model  as  evaluated  by  Kluge/Cain 
Integration  to  find  L  for  a  given  geographic  point.  This  code  will  simply  be  referred  to  as 
Interpolation  to  avoid  confusion  with  the  other  Kluge/Cain  package. 
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The  models  were  compared  primarily  to  determine  if  there  are  experimentally  significant 
differences  which  may  affect  how  the  data  should  be  handled  from  CRRES.  The  primary 
model  currently  in  use  is  the  International  Geomagnetic  R^erence  Field  Revision  1985 
(IGRF  1985).  However,  a  more  recent  paper  by  Cain  reports  a  "jerk"  (discontinuity)  in  the 
secular  variation  of  the  field  in  1983.  This  is  not  accounted  for  in  the  secular  variation 
coefficients  of  the  IGRF  1985  model.  Thus,  Cain  has  ?  wed  a  model  for  epoch  1990  which 
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does  account  for  this  jerk.  In  the  magnetic  field  and  L-shell  comparisons,  no  significant 
effect  was  seen  due  to  this  jerk-  Cain’s  model  agrees  quite  well  with  IGRF 1985  extrapolated 
to  1990.  The  discrepancy  arising  from  this  jerk  should  become  more  apparent  in  future 
epochs.  However,  by  the  time  this  occurs,  IGRF  will  have  been  revised  with  data  extending 
beyond  that  which  Cain  used.  Thus,  this  question  will  be  resolved.  For  the  time  period 
considered  here,  the  difference  is  not  found  'o  be  significant. 

'Hie  software  was  evaluated  to  determine  the  most  efficient  code.  The  GL  codes  use  various 
numerical  techniques  to  integrate  over  the  field  line  in  Cartesian  coordinates  to  determine  B 
and  L  These  techniques  have  been  reviewed  for  accuracy  and  speed  for  both  the  SCATHA 
and  CRRES  satellite  projects.  The  Kluge/Cain  work  uses  inverse  coordinates  centered  on 
the  dipole  so  that  the  integration  is  performed  over  a  straight  line  rather  than  a  curve.  This 
significantly  reduces  the  computation  time.  Interpolation  reduces  it  even  further.  However, 
the  Interpolation  routine  as  it  now  stands  is  only  valid  for  the  epoch  for  which  the  tables 
were  generated.  In  order  to  obtain  the  same  flexibility  in  epoch  that  one  has  with  the 
integration  techniques,  tables  for  other  epochs  must  be  generated  and  a  method  of 
interpolation  between  the  tables  must  be  used. 

The  Interpolation  routine  was  found  to  be  the  fastest.  However,  the  tables  were  determined 
using  Klugc/Cain  Integration  which  incorporated  a  twenty  year  old  value  of  the  dipole 
moment.  The  dipole  moment  is  decreasing  over  time.  In  the  course  of  two  decades,  it  has 
dropped  enough  to  produce  significant  differences  in  the  L-sheil  values.  Thus  for  1990,  it 
would  be  useful  to  generate  interpolation  tables  using  the  current  dipole  moment.  Then  one 
would  be  able  to  fake  advantage  of  this  rapid  technique  to  find  L  for  a  large  set  of  data 
points. 

The  results  from  Kluge/Cain  Integration  agree  well  with  those  from  OPTRACE  and  HMIN. 
It  is  a  faster  routine  than  those  currently  in  use  at  GL  It  is  possible  to  use  this  technique 
with  external  model  routines.  So  for  cases  where  only  a  few  points  arc  required,  or  when 
one  is  interested  in  the  magnetic  field  values,  or  until  the  revised  interpolation  tables  are 
produced,  this  may  be  the  code  to  use. 
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DESCRIPTIONS  OK  THE  MODELS 


2.1  OVERVIEW 

Ibe  four  models  which  are  compared  are  IGRF  1985,  G*in  1990,  MAGSAT  1980,  and 
Burraclough  1975,  IGRF  1985  is  used  ns  the  standard  model  to  which  the  others  arc 
compared,  'litis  is  because  it  is  the  model  released  by  the  International  Association  of 
Geomagnetism  ami  Acronomy  and  is  the  one  which  is  most  widely  used  tor  current  epoch*. 
Ciin  1990  is  considered  here  because  of  a  possible  discontinuity  its  authors  found  in  the 
MMilur  variation  in  1983.  ’Ibis  discontinuity  was  not  incorporated  in  the  IGRF  1985  ntodcl. 
Hence,  as  one  approaches  1990,  significant  errors  may  arise  from  IGRF  1985  due  to  its 
handling  of  the  secular  variation.  MAGSAT  1980  is  presented  because  the  basis  of  this 
model  is  satellite  data.  'Hie  MAGSAT  satellite  was  specifically  designed  for  tltc  pur|x>sc  of 
gathering  high  quality  data  with  good  gtobal  coverage  to  be  used  in  determining  (lie  Barth's 
magnetic  field.  Previous  models  have  all  been  developed  using  ground-bused  data.  Note,  the 
OGO  series  of  satellites  in  the  1960‘s  were  used  in  modeling  the  Sckl,  but  not  to  the  same 
extent  as  MAGSAT.  Burraclough  1975  is  presented  because  it  is  n  model  which  has  been 
used  at  GL  for  years.  '11ms,  it  Is  of  interest  to  those  who  have  used  it  in  the  post  ami  are 
familiar  with  it. 

IGRF  1985,  Giin  1990,  MAGSAT  1980,  ami  Barmc tough  1975  are  all  models  determined  l»y 
tilting  datu  (ground-based  ami  satellite)  to  a  mathematical  model.  Gauss  (1839)  found  that 
the  Barth's  magnetic  field  could  be  expressed  as  the  negative  gradient  of  a  harmonic  scalar 
potential.  This  representation  of  the  field  will  be  discussed.  ’Ilicn  the  data  used  to  obtain 
the  four  models  will  be  briefly  dcscrilted. 

Note,  there  are  other  ways  to  represent  the  field.  Stem  (1976)  reviews  five  groups  of  field 
representations:  a)  those  based  on  current  density,  b)  those  using  a  magnetic  scalar  potential 
(c.g.  Gauss),  c)  poloidal  and  toroidal  fiekfc,  d)  Euler  potentials,  and  e)  local  expansions  of 
the  field  given  a  reference  (mint.  In  a  subsequent  review  paper  by  Backus  (1986),  a  strong 
case  is  made  for  using  poloidal  ami  toroidal  fields  (the  Mic  representation)  rather  than  (tie 
Gaussian  representation.  This  is  because  Gauss  assumes  a  current-free  field.  However,  wlien 
using  satellite  data,  this  assumption  is  not  valid.  lijima  and  Potemra  1 1976)  observed  fieki 
aligned  currents  at  800  km  in  tltc  auroral  zone.  Even  on  quiet  days,  these  led  to  magnetic 
perturbations  of  several  hundred  nanotesla.  Backus  presents  in  detail  the  Mie  representation 
of  the  magnetic  field  in  a  spherical  shell.  'Ibis  particular  representation  lends  itself  to 
combining  satellite  data  with  ground-bused  observations. 

However,  the  Gaussian  representation  is  still  the  most  widely  used  tor  internal  field  models. 
For  low  latitudes  and  altitudes,  the  Gaussian  representation  of  the  internal  fields  b  adequate. 
’Ibe  four  models  considered  in  thb  report  are  all  spherical  harmonic  expansions  of  the 
magnetic  scalar  potential.  Ibe  alternate  representations  mentioned  may  be  of  interest, 
particularly  when  considering  external  field  contributions  as  well  as  the  internal  ones,  but  are 
not  the  focus  of  thb  re|>ort. 
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2.1.1  Derivation  of  the  Main  Field  and  Its  Components 

Assume  the  main  field  can  be  expressed  as  the  negative  gradient  of  a  scalar  potential  (V) 
which  is  a  solution  to  Laplace’s  equation  V2V  -  0.  Then  in  spherical  polar  coordinates,  this 
solution  may  be  written  as  a  series  expansion  in  spherical  harmonics: 

N  n 

V  -  a  £  (a/r)"+1  2  (85  ooa  mX  +  h$J  sin  mX)  I*J(eos  0) 
n-1  m-0 

where 

a  a*  mean  Earth  radius 

r  «  distance  from  Earth's  center 

6  «  geocentric  colatitude  (measured  from  the  geographic  north  pole) 

X  *  longitude  (measured  east  from  Greenwich) 

I**  *  Legendre  polynomial  of  degree  n  and  order  m 

gj,  h*  *  spherical  harmonic  coefficients  that  constitute  the  desired  models 

Note,  several  different  normalizations  are  used  with  the  Legendre  polynomials.  The  models 
presented  here  use  Schmidt  normalization.  This  correlates  the  magnitude  of  a  term  with  its 
contribution  to  the  field.  Hence,  the  coefficients  decrease  with  increasing  n  since  the  lower 
order  terms  contribute  the  most  to  the  field  (at  least  near  the  Earth’s  surface).  There  is  no 
systematic  variation  with  m  for  any  given  n.  Other  normalizations  may  be  obtained  by 
multiplying  the  Legendre  polynomial  by  some  factor,  then  dividing  the  corresponding 
coefficient  by  the  same  factor.  Thus,  V  remains  unchanged  [Stern,  1976]. 

One  may  then  get  the  geocentric  components  of  the  main  field,  B--VV,  by  partial 
differentiation  of  V: 


X  -  Bfl-  (l/r)(5V/B8) 

Y  -  B^-  (-1/r  sin0)(6V/5X) 

Z  -  -Bf-  8V/5r 

where,  X  is  the  northward  component,  Y  is  the  eastward  component,  and  Z  is  the  radially 
inward  component  of  the  field.  The  other  geomagnetic  components  are  related  to  X,  Y,  and 
Z  by: 


D  «■  declination  -  arctan(Y/X) 

I  -  inclination  «  arctan(Z/H) 

H  «  horizontal  intensity  ■»  (X2  +  Y2),/2 
F  «  total  intensity  -  (X2  +  Y2  +  Z2),/2 

In  addition  to  finding  the  main  field  at  a  given  time  (epoch),  one  is  also  interested  in 
extrapolating  any  given  model  to  other  epochs.  To  do  this,  a  secular  variation  (SV)  model  is 
required,  'lliis  is  found  by  assuming  the  coefficients  have  a  linear  dependence  on  time,  eg.  a 
«  a0  +  at,  where  a0  is  the  coefficient  at  some  initial  time  and  t  is  in  years.  In  some  cases,  a 
quadratic  term  (secular  acceleration)  is  also  included  [Stern,  1976]. 
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2.1.2  Model  Epochs  ard  Degrees  of  Derivation 

Variations  between  the  models  are  due  to  two  factors:  a)  the  data  used,  and  b)  the  number  of 
coefficients  determined.  Each  of  these  models  is  centered  on  a  different  epoch,  Barraclough 
1975,  MAGSAT  1960,  IGRF  1965,  and  Cain  1990.  They  have  also  been  derived  to  different 
degrees,  Barraclough  N-12,  MAGSAT  N-66,  IGRF  N-10,  and  Cain  N-15. 

"ihe  distinguishing  qualities  of  the  data  are  when  it  was  taken  and  how  it  was  taken.  'Che 
epoclis  given  with  each  model  reflect  the  date  the  model  is  best  suited  for.  For  MAGSAT, 
the  data  used  was  centered  around  1960.  For  Barraclough,  the  data  was  centered  around 
1%5  then  the  model  was  projected  to  1975.  Similarly,  IGRF  and  Cain  were  based  on  data 
spanning  the  late  1970s  to  the  early  1980s.  ‘Hrese  were  then  projected  to  1985  and  1990. 
Each  of  the  models  uses  a  combination  of  data  sources  including,  permanent  observatories, 
repeat  stations,  ship-towed  magnetometers,  and  satellite  data.  Satellite  data  has  only  been 
incorporated  for  the  past  two  or  three  decades,  lire  satellite  MdataM  included  in  the 
Barraclough  model  was  actually  generated  from  another  model  which  was  based  on  data 
from  the  OGO  scries  of  satellites.  The  quality  of  the  data  varies  due  to  the  global  coverage 
and  the  instrumentation  used.  For  the  ground-based  data  sets,  the  permanent  observatory 
data  is  more  accurate  than  the  ship  data.  However,  in  order  to  adequately  cover  the  globe, 
ship  data  is  essential.  Most  of  the  models  correct  for  the  variation  in  quality  by  using 
weighting  factors.  'Hie  MAGSAT  satellite  data  is  of  particularly  high  quality.  ’Ibis  satellite 
was  designed  specifically  to  measure  the  magnetic  field.  The  one  drawback  to  this  data  is 
that  it  only  covers  a  time  span  of  about  six  months.  Thus,  one  must  be  careful  that  temporal 
variations  in  the  data  don’t  appear  to  be  spatial.  These  data  are  supplemented  with  ground- 
based  data  specifically  to  address  this  problem. 

‘Hie  other  factor  contributing  to  differences  between  the  models  is  the  degree  and  order  to 
which  the  models  are  derived.  The  higher  the  degree  of  the  model,  the  more  complex  the 
structure  becomes.  The  core  component  of  the  field  is  seen  primarily  in  the  lower  degree 
terms  and  is  the  dominant  component  to  about  n-13.  The  higher  degree  terms  reflect  the 
crustal  component  to  the  field.  These  become  dominant  around  n-15.  Thus,  to  most 
accurately  model  the  internal  field  with  all  its  local  variations,  a  model  which  includes  high 
degree  terms  is  desirable.  However,  it  may  also  be  desirable  to  isolate  the  core  component 
of  the  field  and  neglect  smaller  features  due  to  the  crust.  In  this  case,  one  would  truncate 
the  coefficients  at  n£l3.  IGRF  1985  is  truncated  at  n-10  for  this  reason.  Barraclough  and 
Giin  arc  limited  to  their  respective  degrees  (12  and  15)  due  to  the  global  coverage  and 
accuracy  of  iheir  data.  MAGSAT  has  such  a  high  degree  due  to  the  superior  global  coverage 
of  the  satellite,  combined  with  highly  accurate  observatory  data  and  less  accurate  scalar  data 
available  from  remote  regions. 


2.1.3  High  Latitude  Limitations  of  the  Models 

Effort  has  been  made  to  extend  these  models  to  high  latitudes,  but  there  are  problems  in 
doing  so.  The  models  are  derived  by  fitting  the  data  to  the  mathematical  model.  Thus,  to 
obtain  a  reasonable  model  at  high  latitudes,  one  needs  high  latitude  data.  Clearly,  data 
coverage  in  these  regions  is  not  as  complete  as  in  lower  latitude  regions.  IGRF  used  vector 
data  to  ±50°  and  scalar  data  poleward  of  this.  MAGSAT  and  Cain  used  data  to  within  7°  of 
the  poles,  then  interpolated  over  these  regions.  Barraclough  used  data  to  ±85°.  Thus,  one 
must  be  careful  when  applying  these  models  at  high  latitudes. 
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2.2  IGRF  1965 


llte  International  Geomagnetic  Reference  Field  Revision  1985  (hereafter  1GRF  1985)  is  the 
product  of  the  International  Association  of  Geomagnetism  and  Aeronomy  (IAGA)  Division 
I,  Working  Group  1.  This  13  member  group  has  produced  a  series  of  mathematical  models 
of  the  Earth's  main  geomagnetic  field  and  its  secular  variation.  First  adopted  in  1968  and 
now  with  the  fourth  revision  there  nre  five  1GRF  models  (1945,  1950, 1955,  1960  and  1965), 
four  "definitive"  models  (DGRF  1965,  1970,  1975,  and  1980),  and  a  predictive  model 
projecting  the  secular  variation  1985  to  1990  \IAGA,  1986).  The  models  which  will  be 
discussed  here  are  DGRF  1980,  the  predictive  SV  model  1985-1990,  and  IGRF  1985.  Note 
that  the  IGRF  1985  model  was  produced  by  applying  five  years  of  secular  variation  (1980- 
1985)  to  DGRF  1960  [Barraclov.gk,  1967). 

In  the  fourth  revision  of  IGRF,  23  models  from  five  groups  were  submitted  for  models 
Ringing  from  1945  to  1990.  For  the  models  of  interest  here,  the  Working  Group  truncated 
the  model  GSFCMF80  (Goddard  Space  Flight  Center  Main  Field  1960)  by  Langel  and  Estes 
(1985)  to  n-10  and  kept  it  as  the  DGRF  i960.  The  secular  variation  model  for  19S0-19S5 
was  the  mean  of  models  BGSSV82  (British  Geological  Survey  Secular  Variation  19S2)  by 
Quinn  et  al.  [1987],  IZMSV82  (Institute  of  Terrestrial  Magnetism,  Ionosphere  and  Radio 
Wave  Propagation  (IZMIRAN)  Secular  Variation  1982)  by  Golovkov  and  Kolomiitseva 
(1987),  and  U3GSSV87  (United  States  Geological  Survey  Secular  Variation  1967)  by  Peddle 
nnd  Zunde  (1987a).  The  predictive  model  for  1985-1990  is  the  mean  of  BGSSV87 
[llarraclough  and  Kcrridge.  1987),  IZMSV87  [Golovkov  and  Kolomiilteva,  1987),  and 
USGSSV87.  Since  the  IGRF  1985  coefficients  were  determined  by  taking  DGRF  1980  and 
applying  the  secular  variation  for  1980-1965  to  its  coefficients,  the  derivation  of  the  Langel 
and  Estes  model  for  1980  will  be  described. 


2.2.1  'Ihe  Main  Field  Model 

Langel  and  Estes  used  the  data  from  MAGSAT  (November  1979  to  April  1980)  plus  data 
from  91  magnetic  observatories  (1978  to  1982)  for  their  model  of  constant  coefficients 
(degree  and  order  13)  and  secular  variation  coefficients  (degree  and  order  10). 

MAGSAT  had  an  altitude  range  of  300-500  km  with  a  sun-synchronous  orbit  at  the 
terminator  (i.e.  day-night  boundaty).  Thus,  data  was  only  taken  at  dawn  and  dusk  in  local 
time,  'litis  was  the  first  global  vector  survey  of  the  near-Earth  geomagnetic  field.  However, 
in  modeling  the,  field  using  the  gradient  of  the  scalar  potential,  the  field  is  assumed  to  1“ 
curl-free.  This  does  not  hold  at  satellite  altitudes.  MAGSAT  passes  through  field-aligned 
currents  which  have  an  effect  on  the  field  magnitude.  To  compensate  for  this,  component 
data  was  used  only  between  ±50°  geomagnetic  latitude.  Above  this,  only  magnitude  data  was 
used. 

Langel  and  Estes  also  investigated  local  time  asymmetries  by  modeling  the  dawn  and  dusk 
data  separately.  They  found  significant  differences  between  these  two  sets  which  they 
attributed  to  an  eastward  equatorial  electrojet  which  lies  below  the  satellite.  The  electrojet 
affected  the  dusk  data  at  certain  latitudes,  but  not  the  dawn  data,  Thus,  the  dusk  data  was 
corrected  for  B  and  X  between  ±20°  geomagnetic  latitude,  Y  between  ±15°,  and  Z  between 
±50°  using  the  difference  between  the  dusk  and  dawn  potentials.  This  removes  some  of  the 
independence  of  the  dusk  and  dawn  data.  A  model  of  the  electrojet  itself  would  have  been 
the  best  way  to  correct  the  data,  but  one  was  not  available. 
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Hie  cbm  were  put  mio  three  temporal  groups:  November* December  1979,  January-Fcbruary 
I9.v0.  and  March-April  1960.  'Hie  dawn  and  dusk  data  were  separated  and  put  into  5d  x  5° 
equiangular  bins.  Residuals  were  computed  from  a  previous  rr,i*dd  and  data  which  fell 
outside  the  acceptable  range  from  these  were  rejected.  Data  wv  selected  from  each  bin 
such  that  there  were  an  equal  number  of  points  for  equal  areas  all  over  the  globe. 
Preference  was  given  to  dawn  data  and  uncorrected  dusk  data  w'iire  fusible,  with  the 
corrected  dusk  data  used  ns  necessary.  Data  were  also  selected  such  that  a  good  spread  of 
Dst  values  between  ±22.5  nT  was  obtained  to  improve  tire  analysis  of  the  external  rleld.  The 
data  were  estimated  to  be  accurate  to  6nT  root-sum-squarc  (m)  in  each  component  s  id  to 
2riT  rss  in  the  field  magnitude. 

In  addition,  observatory  dam  from  91  stations  (1978-1982)  was  obtained  from  NOAA 
National  Geophysical  and  Solar  Terrestrial  Dam  Center  to  supplement  the  MAGSAT  (him. 
lire  data  were  convened  to  X,  Y,  and  Z  in  geodetic  coordinates  using  an  equatorial  radius  of 
C»mi6S  km  and  flattening  of  298.25. 

For  each  coefficient  the  standard  error  was  determined  as  was  the  error  between  dawn  and 
corrected  dusk  models.  As  in  any  large  statistical  study,  two  sources  of  error  are  present,  one 
in  the  dam  and  the  other  in  the  averaging  process.  The  errors  derived  only  from  the  fit  tend 
to  underestimate  inaccuracies  due  to  model  validity,  improper  estimates  of  dam  correlation, 
and  systematic  effects  of  non-core  fields.  Since  the  dawn  set  and  dusk  set  have  identical 
error  characteristics  except  for  the  effects  of  the  external  fields,  the  difference  in  their  errors 
reflect  external  influences  while  the  coefficient  standard  errors  reflect  internal  inconsistency. 


2.2.2  Tlte  19S0-1985  Secular  Variation  Model 

A  secular  variation  model  was  then  applied  to  DGRF  to  obtain  IGRF  1985.  Three 
candidate  models  were  submitted  to  the  IAGA  committee.  Rather  titan  selecting  one  model 
over  another,  the  committee  decided  to  take  the  mean  of  all  three  and  adopt  that  os  the 
1980-1985  secular  variation  model.  The  candidate  models  are  BGSSV82,  IZMSV82,  and 
USGSSV87.  Each  of  these  will  be  briefly  summarized  below. 

"Hie  BGSSV82  model  [tjainn  el  al.,  1987j  was  based  on  magnetic  observatory  annual  means 
primarily  from  the  World  Digital  Data  Centre  Cl  in  Edinburgh,  supplemented  with  dam 
from  World  Data  Centres  A  and  B.  Annua!  means  from  172  observatories  over  the  period 
1969.5  to  1983.0  were  used.  The  secular  variation  was  estimated  by  finding  the  first 
differences  of  the  component  data.  Linear  fits  to  the  data  were  computed  using  least  squares 
with  weights  varying  linearly  from  1969.5  (given  a  weight  of  0.5)  to  1983.0  (given  a  weight  of 
1.0).  They  found  that  a  single  line  was  inadequate  to  fit  the  data.  This  is  due  to  a 
discontinuity  (or  “"jerk”)  in  tlte  secular  variation  in  1978.  Thus,  they  used  two  line  segments 
to  fit  the  data.  They  obtained  a  model  to  degree  n-8  with  rms  residuals  of  5.8 
nanoTesla/annum  in  X,  4.9  nT  a*1  in  Y,  and  5.2  nT  a*1  in  Z. 

The  IZMSV82  model  [Golovkov  and  Kolomiilteva,  1967]  was  based  on  data  from  160 
observatories.  To  this  they  added  data  at  34  locations  (generally  ocean  areas)  to  fill  gaps  in 
the  data  due  to  the  distribution  of  these  observatories.  Ibe  added  data  was  generated  from  a 
Ijrevious,  model,  1ZMSVSQ.  They  graphically  smoothed  their  X,  Y,  and  Z  data  to  obtain  X, 
Y,  and  Z.  Using  this  technique,  they  tried  to  account  for  some  of  the  non-Iinearitjes  in  the 
data.  They  found  rms  residuals  of  4.0  nT  a4,  3.5  nT  a4,  and  4.6  nT  a4  for  X,  Y,  and  Z, 
respectively. 
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The  USGSSV87  mode!  [Peddie  and  Zunde,  1987*)  was  based  on  data  from  148 
observatories  over  the  period  1980  to  1962.  The  data  consisted  primarily  of  annual  all-day 
means  (i.e.  the  mean  of  all  the  daily  means  for  the  year).  When  it  was  available  for  the 
complete  scries,  quiet-day  means  were  used.  These  are  the  annual  means  based  on  the  five 
international  quiet  days  per  month.  'ihese  are  not  available  from  most  observatories.  It  is 
preferable  to  use  the  quiet-day  means  to  reduce  the  external  contributions  to  the  secular 
variation  models.  A  linear  fit  was  made  to  the  X,  Y,  acd  Z  components,  using  unweighted 
least  squares.  'Ihe  slopes  of  these  lines  were  taken  as  the  estimates  of  X,  Y,  and  Z.  The 
spherical  itarmonic  analysis  was  then  performed  on  these  values.  A  model  of  degree  nn8 
was.  found  to  be  a  gpod  compromise  between  model  size  and  goodness  of  fit  with  6.1  nT  a’1 
in  X,  5.3  nT  a’1  in  Y,  and  5.2  nT  a’1  in  t- 


12.3  The  1985-1990  Secular  Variation  Model 

The  secular  variation  model  adopted  by  the  IAGA  committee  was  also  a  mean  of  three 
models.  'Ihese  three  were  essentially  by  the  same  authors  as  the  models  described  previously 
for  the  1980-1985  SV  model.  'Ihese  models  are  BGSSV87,  IZMSV87,  and  USGSSVS7.  The 
USGSSV87  model  was  discussed  in  the  previous  section.  The  other  two  models  will  be 
briefly  described  below. 

'ihc  BGSSV87  model  [Darraclougk  and  Kerridge,  1967]  was  obtained  from  data  from  159 
observatories  taken  from  1969.5  to  1983.0  .  A  linear  fit  to  the  data  was  made  using  weighted 
least  squares.  'IVo  line  segments  were  required  to  fit  the  data  due  «o  the  jerk  in  1978.  .This 
was  the  same  procedure  as  that  used  to  generate  BGSSV8Z  'Ihen  they  synthesized  X,  Y, 
and  Z.  for  1967  by  extrapolating  the  linear  fit  assuming  a  constant  secular  acceleration  from 
1969.5  to  1987.5.  "Ihis  extended  data  set  was  then  used  to.  ge.t  a  model  to  degree  n-8  with 
rms  residuals  of  10.4  nTa*1, 10.1  nTa*1,  and  9.6  nT  a-1  in  X,  Y,  and  Z,  respectively 

'Ihe  IZMSV87  model  [Golovkov  and  KolomiUtcva,  1987]  used  the  IZMSV82  data  set  from 
160  observatories..  Then,  the  predictive  model  was  obtained  by  using  a  graphical 
extrapolation  of  X,  Y,  and  Z.  They  again  filled  the  data  gaps  in  the  ocean  regions,  but  this 
time  they  used  their  1982  SV  model  to  generate  the  data.  'Ihis  was  to  obtain  a  1985 
predictive  model,  this  in  turn  was  used  to  generate  the  data  to  fill  in  the  gaps  to  obtain  the 
1987  predictive  model. 


2.2.4  Discussion 

As  pointed  out  by  Golovkov  and  Kolomiitseva  [1987],  the  SV  models  for  1985-1990  all  have 
very  similar  background  data.  'Ihus,  differences  in  the  models  arise  due  to  how  the  non- 
uniform  distribution  of  the  observatories  is  handled  and  how  the  data  is  processed.  'Ihese 
differences  can  be  significant  and  can  be  indicative  of  the  reliability  of  the  model.  These 
authors  compared  the  three  candidate  models  for  1985-1990  and  found  that  I21MSV87  and 
USGSSV87  agreed  quite  well,  whereas  the  BGSSV87  model  was  not  in  such  good  agreement. 

Similarly,  Peddie  and  Zunde  [1987b]  looked  at  the  three  candidate  SV  models  for  1985-1990 
and  compared  them  to  data  from  148  magnetic  observatories.  They  assume  that  the  present 
rates  of  change  will  continue.  Ihus,  their  comparison  was  based  on  the  most  recent  trend  in 
SV  (1980-1983  at  these  observatories).  They  also  found  IZMSV87  and  USGSSV87  agreed 
fairly  closely  to  the  data  whereas  BGSSV87  did  not  agree  as  well. 
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Clin  and  Kltith  [1987]  also  reviewed,  these  models.  They  looked  at  the  differences  of  the 
coefficients  and  the  components  X,  Y,  and  Z  between  these  Again,  BGSSV87  was 

fo'ind  to  have  poorer  agreement  with  the  other  two  models  than  they  have  with  each  other. 
BGSSV87  differed  from  the  other  two  by  up  to  100  nT  a*1,  whereas  IZMSV87  and 
USGSSV87  agree  to  within  20  nT  sr1  everywhere  except  in  the  eastern  Pacific 

Rather  than  selecting  one  model  over  another  for  the  SV  model,  the  three  candidate  models 
were  averaged  together  to  comprise  the  SV  1985-1990  model.  The  three  models  were 
apparently  combined  equally  without  any  statistical  weights.  It  is  not  clear  why  this  was 
done.  No  explanation  was  given  by  Barradough  [1987]  in  the  report  which  describes  what 
the  IAGA  committee  selected  (note,  Barradough  was  the  chairperson  of  this  committee). 


2.3  CAIN 

While  conducting  an  analysis  of  candidate  models  for  IGRF  1985,  Cain  and  Kluth  [1987] 
found  what  appeared  to  be  a  jerk  in  the  secular  cliange  determined  from  the  observatory 
data  in  mid-1983.  If  (his  jerk  is  real,  a  question  arises  as  to  how  well  the  IGRF  1985  model 
might  project  to  1990.  They  concluded  that  the  observatories  gave  inadequate  coverage  to 
ascertain  whether  this  was  a  real  phenomenon  or  not.  In  fact,  they  questioned  whether  the 
SV  model  based  on  these  data  alone  give  a  reliable  global  estimate  of  SV  at  all.  Thus,  they 
decided  to  construct  anotlicr  secular  variation  model  based  on  the  observatory  data  plus, 
scalar  data  from  ship-towed  magnetometers  and  data  from  Project  Magnet  aircraft  (from  the 
U.  S.  Naval  Oceanographic  Office). 

They  investigated  linear  and  parabolic  secular  change  models.  As  a  starting  point,  they 
truncated  their  n-66  MAGSAT  model  [Cain  et  a/.,  1968]  to  n-15.  This  was  done  due  to 
the  limitations  of  the  observatory,  Project  Magnet,  and  ship-towed  magnetometer  data  sets. 
These  do  not  have  adequate  global  coverage  to  generate  a  reliable  model  of  such  high 
degree.  These  additional  data  sets  were  used  to  determine  a  secular  change  model  to 
extrapolate  MAGSAT  to  future  epochs.  Various  secular  variation  models  were  derived  with 
a  linear  model  (M8386L)  giving  the  best  projection  of  the  field  to  1990  (Cain,  1987]. 

Note,  in  this  section,  only  the  derivation  of  the  secular  change  model  which  is  applied  to  the 
MAGSAT  model  to  obtain  the  Grin  1990  model  is  discussed.  The  derivation  of  the 
MAGSAT  model  which  is  used  here  as  the  main  field  model  will  be  presented  in  the 
subsequent  section  of  this  report  (2.3  MAGSAT). 


2.3.1  Description  of  the  Secular  Variation  Data 

The  observatory  data  came  from  several  sources,  digital  data  at  one  and  2.5  minute  intervals 
and  tabular  data  of  hourly  and  annual  means  with  absolute  scale  offsets  (most  of  the  time). 
Most  of  the  data  was  supplied  by  the  National  Oceanic  and  Atmospheric  Administration 
(NOAA)  in  Boulder,  Colorado.  Additional  data  was  obtained  from  the  World  Data  Center 
B  in  Moscow.  Plus,  some  data  was  sent  directly  from  the  magnetic  observatories. 

The  observatory  data  needed  to  be  put  into  a  uniform  format.  Thus,  they  used  hourly 
overages  when  it  was  possible  and  supplemented  these  with  the  annual  all-day  means.  Cain 
[1987]  points  out  that  the  quiet-day  means  are  preferable  to  the  all-day  means  (the  former 
determine  the  annual  mean  using  only  the  means  from  the  five  international  quiet  days, 
whereas  the  latter  uses  all  of  the  daily  means  to  evaluate  the  annual  mean).  However,  these 
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are  not  as  readily  available  os  the  annual  means  determined  from  the  all-day  means.  Data 
was  selected  according  to  following  criteria.  To  maximize  the  amount  of  data  while 
minimizing  Kp,  Kp<2  was  chosen  as  the  cutoff  for  acceptable  data.  Similarly,  to  reduce  the 
effects  of  diurnal  variations,  only  data  within  four  houn  of  local  midnight  was  kept.  To 
account  for  external  contributions  to  the  field,  the  Dst  correction  was  calculated  from  near 
equatorial  H  values.  There  has  not  been  a  global  determination  of  the  zero-level  Dst;  thus 
they  assumed  it  is  constant  over  the  data  set  and  that  its  average  may  be  derived.  To 
simplify  things  somewhat,  a  ratio  of  internal  to  external  effects  equal  to  0.28  was  used,  ‘litis 
may  be  too  high  for  Kp<2  and  should  be  reconsidered.  All  of  the  data  was  expressed  in 
terms  of  the  X,  Y,  Z  components,  with  the  incomplete  vector  data  dropped,  and  each 
component  then  corrected  for  Dst. 

‘Hie  Project  Magnet  data  were  provided  by  NOAA  and  by  the  Navy  in  five  minute  intervals. 
It  was  binned  by  year  and  by  areas  equivalent  to  an  equatorial  10°  x  10°  square.  Only  data 
with  Kp^2  were  retained.  Data  within  two  hours  of  local  noon  were  deleted.  Comparisons 
were  made  to  a  previous  model,  with  data  differing  from  the  model  by  more  than  400  nT 
rejected.  Vector  data  was  weighted  and  used  along  with  the  more  accurate  scalar  data.  Dst 
was  corrected  for  in  the  same  manner  ns  the  observatory  data. 

‘ilic  ship-towed  data  is  a  vast  database.  Hence,  one  could  be  selective  about  which  data  to 
retain.  Values  were  selected  for  Kp<l+  and  within  four  hours  of  local  midnight.  Again  the 
Dst  correction  was  the  same  as  the  rest  and  the  binning  was  by  year. 

All  the  data  was  combined  into  a  40.000  record  set.  All  the  observatory  data  was  used  with 
the  Project  Magnet  (vector  and  scalar)  and  ship  towed  (scalar)  data  filling  in  the  rest. 


2.3.2  Derivation  of  the  Secular  Variation  Model 

To  obtain  secular  variation,  the  X,  Y,  and  Z  components  are  fit  to  a  polynomial  using  least 
squares.  Then  the  slope  at  a  given  time  is  determined  to  obtain  X,  Y,  and  Z.  To  make  the 
calculations  easier,  the  hourly  averages  were  grouped  into  10  day  means.  The  data  set  was 
put  into  several  test  models  to  optimize  weighting  factors  and  to  check  the  stability  of  the 
solutions  with  the  available  data  distributions.  'Hie  test  models  led  to  a  final  parabolic  model 
of  degree  n«7.  To  find  the  best  projection  of  the  field  to  1990,  it  was  decided  to  use  this 
parabolic  model  to  generate  additional  data  to  compensate  for  the  lack  of  available  data 
post-1983.  Thus,  the  real  data  from  i960  to  1983  results  in  the  best  model  for  this  epoch, 
while  the  calculated  "data”  influences  the  secular  variation  after  1983.  It  is  assumed  that  a 
linear  projection  to  future  epochs  is  the  best.  Here,  two  linear  models  were  obtained,  one 
for  1980-1953  (M8083L)  and  the  other  for  1983  to  the  end  of  the  data  set  (M8386L).  'Iliese 
models  were  derived  to  degree  n-7  (stable  solutions  did  not  exist  beyond  this).  Hie  latter 
model  gives  the  best  projection  to  1990. 

Giin  [1987)  states  M8386L  gives  the  best  projection  to  1990.  IGRF  1985  does  not  account 
for  the  mid-1983  jerk,  thus  as  one  gets  further  away  from  1983,  it  is  anticipated  that  Giin 
1990  will  represent  the  field  better  than  IGRF  1985.  Project  Magnet  and  ship  towed  data 
are  consistent  with  the  post-1983  model  and  confirm  the  jerk  detected  in  1983  as  well. 
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2.4  MAGSAT 


MAGSAT  was  a  satellite  which  was  placed  in  a  sun-synchronous  orbit  in  the  dusk-dawn 
meridian  for  the  specific  purpose  of  measuring  the  Earth’s  magnetic  field  (See  Section  2.1.1). 
Several  people  have  used  the  MAGSAT  database  to  derive  magnetic  field  models.  For  this 
report,  coefficients  derived  by  Cain  ct  al.  (1969a)  for  an  n-66  spherical  harmonic  expansion 
have  been  used.  This  high  degree  is  obtained  by  using  the  Neumann  method  of  Gauss* 
Legendre  quadrature  (ScAmits  el  al.,  1989).  Although  the  model  goes  to  n*66,  the  set 
available  at  GL  only  contains  the  coefficients  to  n«50.  For  computational  expediency',  one 
may  select  a  subset  of  this  which  has  been  truncated  to  n-15. 


2.4.1  Description  of  tire  Data 

Gun  cl  al.  (1989a)  applied  the  Neumann  method  of  Gnuss-Lcgcndre  quadrature  to  the  radial 
components  (Bf)  of  the  MAGSAT  data.  Starting  with  a  previous  MAGSAT  model,  new 
observatory  data  was  added.  To  obtain  better  secular  change  (SV)  data,  hourly  values  were 
only  selected  within  two  hours  of  midnight  and  for  a  limited  range  beyond  Kp«l+.  These 
were  then  corrected  for  Dst,  averaged  for  each  universal  time  day  and  converted  to  X,  Y,  2. 
Using  least  squares,  X,  Y,  2  were  found.  The  fit  was  good  except  at  the  North  Pole  due  to 
seasonal  variations  in  the  polar  ionospheric  currents,  which  lead  to  secular  changes  which  are 
higher  than  one  expects  from  the  internal  field. 

In  order  to  minimize  annual  variations  which  arise  in  data  sets  taken  over  a  short  period  of 
time,  data  was  added  to  expand  the  collection  time  of  MAGSAT.  However,  it  was  not 
extended  too  far  so  as  to  avoid  problems  with  the  jerks  in  1978  and  1983.  Data  from 
MAGSAT  was  selected  every  128th  observation  (roughly  50  km  along  satellite  orbit).  Only 
values  where  Kp£2+  and  which  were  within  the  100  nT  of  the  previous  model  were  kept.  Of 
the  1,330,285  observations,  half  were  lost  to  the  Kp  requirement  and  4478  were  lost  to  the 
100  nT  requirement.  Data  were  binned  every  3°  in  longitude  and  into  64  latitudes  (the  roots 
of  the  n-64  Legendre  polynomial).  Using  Sugiura’s  computation  of  Dst,  the  data  near  the 
equator  were  corrected.  Note,  this  does  not  account  for  a  zero-level  in  Dst;  thus,  there  may 
be  an  offset  to  this  correction.  This  computation  was  done  only  for  near  equator  stations 
and  is  not  easily  extrapolated  to  higher  latitudes.  So  beyond  the  range  of  Sugiura’s 
correction,  a  constant  ratio  of  .28  for  the  internal  and  external  terms  was  used  to  correct  the 
data  for  Dst  contributions.  This  was  checked  via  scatter  plots  of  the  correction  versus 
observatory  hourly  values  when  Kp<2+,  and  was  deemed  reasonable. 

In  getting  the  residuals,  the  average  external  term  was  canceled  in  taking  the  differences  since 
the  average  external  term  is  calculated  along  with  the  internal  terms.  However,  the  average 
induced  internal  field  h  still  unaccounted  for.  Five  percent  of  the  data  was  divided  into  dusk 
and  dawn  according  to  the  local  time  of  the  observation  and  averaged  by  dip  latitude. 
Analysis  of  these,  points  to  a  meridional  current  on  the  dusk  meridian  of  the  Y  component 
(A faeda  <l  al.,  1985]  and  to  a  weak  westward  electrojet  flowing  at  dawn  below  MAGSAT 
altitudes  in  the  E-region  [Gain  et  al.,  1989a).  These  were  corrected  for  by  using  a  simple 
function  of  dip  colatitude  to  reduce  the  peak  scatter  (by  as  much  as  15  nT  at  some  latitudes). 


2.4.2  Results  from  the  Model 

Ultimately,  100  to  400  data  values  were  averaged  per  block  with  standard  deviations  of  the 
means  ranging  from  3  nT  to  40  nT  with  most  under  10  nT.  Using  Neumann’s  method  on  a 
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dam  set  with  interpolated  polar  values,  the  same  anomalies  seen  st  n=29  were  still  present, 
but  mere  sharply  focussed  in  this  case.  Beyond  n-35,  some  of  the  features  may  be  of 
geologic  significance  [Cain  el  « l 1963,  Cain  tt  *1.,  1969b],  However,  projections  to  the 
surface  showed  north-south  siriations  which  may  be  the  result  of  inadequate  noise  reduction. 
'Hie  low  degree  terms  in  this  model  were  also  adjusted  to  look  at  the  core  components 
(n»ll  to  M)  which  are  cut-off  in  the  n-10  models  (i.e.,  IGRF  and  DGRF  models).  This 
technique  makes  it  possible  to  obtain  high  degree  models. 

For  the  purposes  of  this  report,  small  sample  fields  for  1960  were  compared  in  see  the 
differences  in  L-$hdl  computed  from  the  n-50  model,  n«50  truncated  to  20,  and  n«50 
truncated  to  15.  The  differences  in  the  tested  areas  were  very  small,  on  the  order  of  10-4  - 
10-5  for  low  altitude  and  no  difference  for  higher  altitude  (>1000  km).  One  sample  point 
showed  a  difference  of  .1-1  in  I.  which  is  significant.  However,  the  increase  in  computational 
time  for  the  n">50  model  did  not  seem  justified  for  these  comparisons.  Thus,  the  model 
truncated  to  n"*15  was  used  in  the  comparisons  with  Barraclough  1975,  IGRF  1985,  anti 
Giin  1990. 


15  BARRACLOUGH 

Barraclough  et  at.  [1975]  derive  three  spherical  harmonic  models,  the  main  field,  the  secular 
variation,  and  the  secular  acceleration.  The  main  field  model  is  composed  of  168  spherical 
harmonic  coefficients  to  degree  and  order  11  The  secular  variation  model  has  80  coefficients 
with  degree  and  order  8.  'Hie  secular  acceleration  model  is  comprised  of  the  26  most 
significant  coefficients. 


2.5.1  General  Description  of  the  Data 

The  main  field  coefficients  are  Schmidt  quasi-normalized  in  nano-Tesla.  They  were  derived 
from  five  data  sets:  survey,  observatory,  oceanographic,  ships*  compass,  and  satellite.  The 
first  four  of  these  include  all  the  1955  to  1975  data  from  the  World  Digital  Data  Centre  Cl 
at  Herstmonceux.  All  the  data  were  reduced  to  sea  level  and  a  common  epoch,  where  the 
common  epoch  was  chosen  to  be  1965.  This  was  selecied  because  it  is  the  mean  date  of  the 
observations  and  thus  as  many  data  are  extrapolated  forward  as  backward,  thereby 
minimizing  SV  errors. 

'Hie  survey  data  were  obtained  from  land,  sea  and  aircraft  magnetic  surveys  of  one  or  more 
of  the  following  field  components:  Declination  (D),  inclination  (I),  horizontal  intensity  (H), 
vertical  intensity  (Z),  and  total  intensity  (F).  The  observatory  data  are  the  annual  mean 
values  from  approximately  200  fixed  stations  with  uneven  global  distribution.  The  observed 
components  are  generally  D,  H,  and  Z  or  I,  but  some  of  the  high  latitude  stations  also 
supplied  north  intensity  (X)  and  east  intensity  (Y).  The  oceanographic  data  are  comprised 
of  a  large  number  of  total  intensity  observations  made  on  surveys  which  covered  most  ocean 
areas.  The  ships’  compass  observations  were  used  to  fill  in  sparsely  covered  regions.  They 
are  generally  D  values  of  limited  accuracy  from  compass  record  books  of  merchant  vessels. 
Finally,  the  satellite  data  were  obtained  from  a  model  (Cain  8/73)  derived  from  the  OGO 
scries  of  satellites  (numbers  2,  4,  and  6,  collectively  referred  to  as  POGO).  The  satellite  data 
were  handled  differently  and  will  be  discussed  later. 

'lire  data  sets  were  binned  into  1654  tesserae  (regions  defined  by  any  two  given  latitude  lines 
and  any  two  given  longitude  lines)  of  approximately  equal  area  (5°  x  5°  at  the  equator,  5°  x 
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120°  at  the  poles).  Unlike  the  other  models,  the  data  for  Barraelough  1975  are  not  centered 
on  1975,  but  rather  on  1955.  Thus,  the  data  were  evaluated  using  1%5.Q  as  the  mean  time  of 
the  data  set,  then  the  model  was  projected  to  1975.0.  The  mean  value  of  the  data  in  each 
tessera  was  found  with  a  separate  series  of  means  for  each  data  set  obtained.  Post-1965  data 
were  reduced  to  19o5.Q  using  a  preliminary  1970  SV  (secular  variation)  model  by  Malm  and 
Dark  (1974).  Uptimes  front  1960  to  1965  were  done  using  the  i960  SV  model  by  Malin 
1 1969)  and  the  1965  SV  model  by  Malin  and  Dark  (1974).  Error  checking  of  the  tcsscr.il 
means  using  residuals  between  them  and  the  corresponding  IGRF  1965  values  were 
performed. 

'Hie  separate  means  for  each  data  set  were  combined  using  a  weighting  system  based  on  the 
number  of  observations,  the  date  of  the  observation  (1955-1965  was  weighted  by  .5  due  to 
their  lesser  importance),  and  the  type  of  data  (compass  data  was  weighted  by  .01  since  it  is 
much  less  accurate),  '(lien  a  differential  method  was  used  [Cain  el  a/.,  1%7)  to  obtain  a 
spherical  harmonic  model  of  the  main  Held  from  alt  of  the  elements  of  the  geomagnetic 
field.  With  this  technique,  one  solves  for  corrections  to  ihc  coefficients  of  an  approximate 
spherical  harmonic  model,  1GRF  1965  in  this  ease. 

At  this  point,  (he  satellite  data  set  was  incorporated  into  the  model.  Rather  than  re¬ 
evaluating  tltc  POGO  data.  Barraelough  et  al.  used  the  TOGO  (8/73)  model  by  Cain.  ‘Ibis 
model  includes  the  main  field  and  secular  variation  to  degree  and  order  14  for  epoch  1970. 
With  this,  total  intensity  values  were  generated  at  5°  intervals  of  latitude  and  longitude 
between  ±85°  at  S00  km  (like  real  satellite  data)  for  epoch  1965.0.  These  2520  values  were 
then  used  to  obtain  normal  equations  similar  to  the  surface  data  equations.  Then  (he 
satellite  set  and  surface  set  were  given  equal  weights  to  obtain  the  spherical  harmonic  model 
based  on  all  (he  data. 


2.5.2  Secular  Variation 

The  secular  variation  model  was  derived  from  annual  means  from  magnetic  observatories, 
repeat  stations  (stations  where  observations  are  less  frequent  than  once  per  day),  survey  data, 
and  satellite  data.  ‘Ihc  magnetic  observatory  data  set  caine  from  180  observatories  and  is 
comprised  of  D,  H,  and  Z  observations.  The  differences  of  each  element  were  plotted  versus 
time  with  the  curves  extrapolated  to  1975.  Each  extrapolated  value  was  given  a  range  in 
which  the  actual  value  might  fall  based  on  (lie  reliability  of  the  data  and  on  correlations  with 
nearby  observatories.  This  uncertainty  range  was  then  used  to  assign  weights  to  the  532 
values.  Again,  extensive  error  checking  was  performed  on  each  data  set  using  the  average  of 
the  other  four  as  the  nprm.  ‘Ihe  best  linear  fit  of  data  from  656  stations  (1877  data  points) 
was  used  to  determine  D,  H,  Z  at  the  nv<m  epoch  of  observation. 

X,  Y,  and  Z  values  were  generated  from  the  difference  between  the  1960  and  the  1970  eighth 
order  main  field  spherical  harmonic  models.  These  were  derived  from  1955-1965  survey  data 
and  post-1965  survey  and  oceanographic  data,  respectively.  These  data  points  were  generated 
at  the  center  of  earu  tessera  and  weighted  according  to  the  number  of  observations  from 
1960  or  1970  in  that  tessera  (the  lower  of  the  two  if  they  are  not  the  same).  Similarly, 
POGO  (8/73)  was  used  to  get  an  interpolating  function  by  synthesizing  F. 

Of  the  four  categories  (observatory,  repeat  stations,  survey  and  satellite  data),  only  the  first 
refers  to  1975,  the  rest  were  updated  with  a  secular  acceleration  model.  A  separate  model 
was  derived  from  each  group  and  the  models  were  compared  to  each  other  to  double  check 
the  suitability  of  the  observational  equations  used. 
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2.5.3  Secular  Acceleration 


‘Hie  secular  acoekration  model  used  seven  SV  models.  Hie  five  year  means  of  X,  Y,  and  Z 
from  the  same  set  of  80  observations  (giver,  equal  weight)  for  four  models  (to  6<h  order  and 
degree)  1942.5-1947.5,  19475-19525,  19525-19575, 19575-19625  were  used  Similarly  data 
from  the  19625-19675  model  was  used,  but  this  model  was  based  on  data  from  118 
observatories.  The  final  two  models,  1970  and  1975,  were  based  on  annual  D,  H,  and  Z 
means  from  180  observatories  weighted  according  to  confidence  limits  on  data.  These  last 
two  models  are  complete  to  5th  order  and  degree,  but  they  alia  have  some  coefficients  out  to 
Sth  degree. 

'11k  linearity  of  the  SA  plots  indicate  higher  order  derivatives  are  negligibk  in  comparison  to 
SA.  The  coefficient*  were  calculated  using  a  least  squares  fit  to  the  coefficients  of  the  seven 
SV  models  of  equal  weight.  Standard  deviations  were  calculated  from  the  scatter  of  points 
about  the  best  straight  line.  The  final  model  then,  is  composed  of  26  coefficients  exceeding 
the  standard  deviation  by  more  than  25. 
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3.  B  FIELD  COMPARISONS  OF  THE  MODELS 


Each  of  the  models  (IGRF  19S5.  Din  1990,  MAOSAT  I960,  and  Barraclough  1975)  were 
used  to  obtain  the  total  B  Field  and  let  component*  at  five  altitudes  (350  km,  850  km,  15(30 
km.  10,000  km,  and  15.000  km)  for  1990,  1985,  and  1900,  The  differences  between  these 
models  were  then  plotted  for  each  altitude  and  year  to  get  an  overview  of  the  model 
discrepancies.  IGRF  1963  »  wed  as  the  sun*'  •!  since  it  is  the  best  known  and  most  widely 
used  of  all  the  models.  More  detailed  comparisons  were  made  by  looking  at  B  values 
obtained  from  each  model  over  a  range  of  latitudes  lor  a  given  longitude,  epoch,  and  altitude. 
Initially,  the  whole  latitude  range  was  used  (±75°).  liven,  the  range  CRRES  will  be  In  was 
examined  more  closely  (±20°).  Finally,  contour  plots  of  the  field  determined  from  IGRF 
1985  for  1960  and  1990  at  350  km  are  presented  to  show  the  overall  change  In  tlie  field  over 
this  ten  year  span. 


3.1  OVERALL  MODEL  DIFFERENCES 

Figures  3.1a  and  3.1b  arc  examples  of  the  contour  plots  produced  initially  to  look  for  trends 
in  die  model  differences  over  Use  full  latitude  and  longitude  ranges  considered.  Figure  3.1a 
shows  the  differences  found  for  1990  at  an  altitude  of  350  km.  *lhe  lower  right  panel  of  the 
plot  shows  the  magnitude  of  B  over  the  globe  as  calculated  by  IGRF  1985.  The  contours  arc 
drawn  at  increments  of  5000  nT.  The  South  Atlantic  Anomaly  (SAA)  is  dearly  visible  as  arc 
the  poles.  Note,  the  appearance  of  the  two  maxima  in  the  northern  hemisphere  are  a 
product  of  the  projection.  If  this  plot  were  mapped  onto  a  sphere,  there  would  only  be  one 
northern  maximum  indicative  of  (he  north  magnets  x>lc.  ‘Hie  magnetic  field  is  strongest  at 
the  poles,  decreasing  towards  the  equator  as  is  exp  -ed  for  a  nearly  dipole  field.  The  other 
three  panels  show  the  absolute  differences  in  nT  between  Din,  MAOSAT,  and  Barradough 
versus  IGRF,  respectively. 

The  differences  seen  in  Din-IGRF  and  MAGSAT-IGRF  are  less  than  one  percent  of  the 
total  field.  There  docs  not  appear  to  be  any  systematic  difference  in  either  of  these  plots. 
'Hie  scattered  distribution  of  the  difference  contours  is  probably  a  reflection  of  the  different 
degrees  of  these  models.  Both  Din  and  MAG5AT  have  higher  order  coefficients  than 
IGRF.  Thus,  they  represent  some  of  the  local  crustal  contributions  to  the  fkkl  which  IGRF 
does  not.  Still  these  are  small  differences.  The  differences  seen  in  Barroclough-IGRF  are 
somewhat  larger  than  those  observed  in  the  other  two  plots.  They  are  still  only  on  the  order 
of  a  percent  or  so.  Here,  the  distribution  of  the  contours  is  not  scattered  as  it  is  with  the 
others.  This  implies  a  more  systematic  variation  from  IGRF  which  is  probably  attributable 
to  the  projection  of  the  model  over  15  years  to  1990. 

Figure  3.1b  is  essentially  the  same  as  Figure  3.1a,  however  it  shows  the  differences  at  an 
altitude  of  10,000  km.  This  altitude  is  about  as  distant  as  one  should  go  without  including 
the  external  field  calculations.  Here,  the  field  strength  has  dropped  off  by  more  than  an 
order  of  magnitude.  The  polar  maxima  are  still  visible,  but  they  are  not  as  distinct.  The 
SAA  is  still  apparent,  but  again  it  is  not  as  distinct  as  it  is  at  perigee. 

The  differences  seen  in  Din-IGRF  and  MAGSAT-IGRF  are  less  than  one  percent.  'Ihe 
contours  are  not  as  scattered  as  previously.  The  crustal  contributions  to  the  field  are  not  as 
pronounced  as  at  lower  altitudes.  Din-IGRF  shows  some  differences  centered  around  the 
SAA  and  some  differences  centered  over  the  southern  Indian  and  Pacific  Oceans.  This  is  to 
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Figure  3.1b.  Contour  plots  of  model  B,.,  differences  in  nanoTesJa  covering  the  full  latitude  and  longitude 
ranges.  IGRF  1935  is  used  as  the  reference  model.  The  comparisons  arc  done  for  an  altitude  of  10,000 
km  at  1990.0. 


be  expected  due  to  the  poorer  data  coverage  in  the  southern  ocean  regions.  Similar  features 
ore  seen  in  MAGSAT-IGRF.  Here,  there  is  no  distinct  difference  centered  on  the  SAA. 
However,  there  are  a  couple  of  distinct  northern  features.  One  set  of  contours  is  centered 
over  the  northern  Pacific  and  another  is  centered  over  western  China  (a  particularly 
mountainous  region,  presumably  this  area  has  not  been  covered  well  with  ground-based 
observations}.  Again,  the  differences  seen  in  both  of  these  plots  are  quite  small  compared  to 
the  total  field  observed.  The  differences  seen  in  Barrackmgh-IGRF  are  larger  than  the 
others,  but  still  less  than  one  percent  of  the  magnitude  of  the  field.  The  largest  differences 
are  centered  around  western  China  and  the  region  around  Baffin  Bay  (northeastern  Canada), 
'litis  latter  area  is  surprising,  since  this  region  is  reasonably  well  covered  by  magnetic 
observatories. 

Table  3.1  presents  the  maximum  differences  between  ±20°  latitude  for  aii  of  the  years  and 
altitudes  calculated.  Each  model  is  compared  to  IGRF  1965  (model-IGRF).  MAGSAT 
shows  the  closest  agreement  at  1980.  By  1990,  Cain  agrees  best  with  IGRF,  For 
Barraclough  and  MAGSAT  the  differences  steadily  increase  with  time.  For  Cain,  the 
differences  are  smallest  for  1965  increasing  as  one  extrapolates  both  forward  and  backward 
from  this  epoch.  Note,  that  while  the  absolute  differences  drop  with  increasing  altitude,  so 
does  the  magnitude  of  the  field.  The  percent  differences  also  drop  with  increasing  altitude 
slightly.  Tlicsc  worst  differences  are  all  less  than  one  percent  of  the  magnitude  of  the  field 
for  1980.  MAGSAT  und  Barraciough  are  off  by  about  1.5%  for  1990  at  350  km.  For  10.000 
km,  these  vary  from  IGRF  by  less  than  a  percent. 


3.2  COMPARISONS  AT  SPECIFIED  LONGITUDES 

Next,  line  plots  of  the  magnitude  of  B  were  generated  to  high  latitude  at  four  longitudes  over 
the  globe  (see  Figure  3.2}.  All  four  models  are  plotted  in  each  panel.  The  agreement  is 
close,  but  there  are  clearly  areas  where  the  models  do  not  coincide.  As  is  expected, 
Bnrraclough  varies  the  most,  but  careful  examination  (particularly  of  ine  lower  left  panel) 
shows  that  MAGSAT,  Cain,  and  IGRF  do  not  completely  overlap  each  other  either.  These 
plots  show  polar  maxima  (the  right  two  panels)  and  the  SAA  minimum  (the  lower  right 
panel).  Note,  that  the  longitudes  shown  are  not  at  the  peaks  of  any  of  these  features,  but 
one  can  see  how  the  field  is  changing  over  longitude.  Taking  a  closer  look  at  this  sort  of 
plot,  Figure  3.3a  focusses  in  on  tne  ±20°  range  at  1990  (upper  panel)  and  1980  (lower  panel) 
for  80°  east  longitude  at  perigee.  The  differences  are  apparent  in  1980,  but  they  have 
significantly  increased  by  1990.  Figure  3.3b  gives  a  similar  comparison,  but  for  an  altitude  of 
10,000  km.  The  models  are  in  much  better  agreement  for  both  epochs  at  this  altitude. 


3.3  B  DETERMINED  BY  IGRF  1985  FOR  1990 

Figure  3.4a  shows  the  difference  between  the  constant  B  contours  calculated  using  IGRF 
19S5  for  epochs  1990  and  1980  at  an  altitude  of  350  km.  While  the  general  shape  and 
distribution  of  the  field  lines  remains  unchanged  over  the  ten  year  span,  the  field  strength  is 
decreasing  (this  is  most  apparent  at  the  edges  of  the  plot).  This  is  seen  to  be  the  case  at  all 
five  altitudes  which  were  compared.  This  is  due  to  the  fact  that  the  dipole  moment  is 
decreasing  (about  0.05%  annually).  The  decreasing  dipoie  has  significant  ramifications  in 
computing  the  magnetic  field  (see  section  4.3.3.1).  Figures  3.4b  and  3.4c  are  included  for  the 
benefit  of  the  reader.  These  show  the  1990  B  magnitude  contours  determined  from  IGRF 
1985  for  the  remaining  four  altitudes:  850  km,  1500  km,  10,000  km,  and  15,000  km. 
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Table  3.1.  MAXIMUM  11^.  DIFFERENCES  BETWEEN  ±20°  LATITUDE 
All  models  Cbmpartd  lo  IGRF  1985  (Model-IGRF) 


Year 

All.  (km) 

Cain  (n‘0 

Magiat  (nT) 

Barraelough  (nT) 

1990 

350 

257.34 

-306.22 

-639.74 

850 

162.7? 

-175.11 

-409.83 

1500 

97.12 

-98.37 

-249.93 

10000 

3.28 

-5.58 

-14.37 

15000 

1.19 

-2.36 

-5.39 

1985 

350 

95.57 

-159.15 

-442.93 

850 

60.16 

-90.13 

-287.49 

1500 

36.33 

-51.14 

-175.70 

10000 

1.59 

-2.65 

-9.50 

15000 

0.61 

-1.13 

-3.54 

1‘>.S0 

350 

163.59 

66.05 

-249.10 

850 

105.47 

40.45 

-165.36 

1500 

63.34 

23.26 

-106.03 

10000 

1.95 

-1.60 

-4.97 

15000 

0.72 

*0.71 

-1.73 
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Figure  3.2.  Line  plots  of  the  four  models  plotted  together  at  four  different  longitudes  over  the  globe  for 
the  full  latitude  range.  For  the  most  part  the  models  agree  at  this  resolution,  although  Barradouglt  is  seen 
to  vaty  somewhat  in  each  of  the  plots. 


Omog  vi.  lotlluda:  Attitude  «  330.00  Epoch  ■  1000.0  Longitude  •*  80 


Bmag  v*.  Latitude:  Altitude  «  330.00  Epoch  ■  1080.0  Longitude  •*  80 


Figure  3.3a.  Line  plots  similar  to  those  in  Figure  3.2,  but  over  a  smaller  latitude  range.  At  this  resolution 
differences  between  the  models  are  apparent.  The  top  panel  shows  the  differences  at  350  km  and  80*  cast 
longitude  for  1990.0.  The  bottom  panel  is  the  same  plot  done  for  epoch  1980.0.  Clearly  discrepancies  arise 
o\cr  time. 


t--r.cs  (nono-Ttiio)  £■■*=$  (oorto-7c»5a) 


LONGITUDE  (DECREES) 

IGRF85  8,^.  Cbnlours  Yeor  ■  1980.0  Altitude  « 


LArnuoc  (oecnccsj  (occnscs) 


Figure  3.4c.  Contours  of  B  as  determined  by  IGRF 1985  for  1990  over  the  full  range  of  latitudes  and 
longitudes  at  10,000  km  (top)  and  15,000  lan  (bottom). 
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DISCUSSION  OF  SOFTWARE 


*1.1  INTRODUCTION 

Three  basic  software  packages  were  used  in  this  analysis.  One  has  been  derived  from 
OPTRACE  [Radcx,  Inc.,  1987],  which  was  developed  to  compute  the  magnetic  field 
quantities  (model  field  vector,  L  parameter,  100  km  field  line  intercepts,  etc)  which  will  be 
stored  in  the  CRRES  ephemeris.  This  package  incorporates  the  IGRF  196S  model,  with  the 
secular  variation,  for  the  internal  portion  of  the  magnetic  field,  and  the  Obon-Pfitter  tilt* 
dependent  model  [OUon  and  PfiUtr,  1977]  for  the  external  portion.  This  package  was 
modified  to  include  lire  other  three  models  studied  in  this  report  (Darraclough  1975,  Magsat 
1980,  and  Giin  1990),  and  to  exclude  the  external  field.  Finally,  to  save  computation  time, 
the  computations  of  parameters  other  than  the  magnetic  field  vector  and  L  were  removed. 
This  final  package,  calleJ  HMIN,  was  tested  to  assure  exact  reproduction  of  OPTRACE  L 
values. 

The  other  wo  packages,  received  from  J.  C  Giin,  are  resurrections  of  software  developed  by 
Kluge  (1972],  and  Kluge  and  Lenhart  (1972),  respectively.  The  first  package,  SHELLC,  uses 
a  very  fast  algorithm  [Kluge,  1970]  for  computation  of  the  L  parameter.  The  second 
package,  INTEL,  interpolates  L  from  a  prestored  table,  and  is  therefore  even  faster. 
However,  the  table  is  valid  only  for  one  epoch,  i.e.,  no  secular  variation  is  included.  The 
table  was  generated  by  Giin  for  his  model  at  the  1990  epoch.  However  we  do  not  at  present 
possess  the  full  software  needed  for  generating  tables  for  other  models  and  epochs.  The 
resurrected  software  was  tested  (Cain,  1987]  with  old  cases  run  by  Kluge  to  assure 
reproduction  of  results. 


<1.2  DEFINITION  OF  THE  L  PARAMEFER 

'Hie  L  parameter  maps  points  in  a  given  magnetic  field  to  points  in  a  reference  dipole  field. 
The  L  parameter  of  a  point  [Mcllwain,  1961]  in  a  given  magnetic  field  is  the  distance  from 
the  dipole,  in  earth  radii,  of  the  equatorial  intersection  of  a  dipole  field  line  passing  through 
a  point  with  the  same  invariant  integral  I  and  field  strength  B  in  the  reference  dipole  field  as 
the  point  of  interest  in  the  given  field.  The  invariant  integral  at  a  point  A  is  the  familiar 
field  line  integral 

I  -  J£(l-B,/Bpds 

where  C  is  the  conjugate  point  to  A,  ds  is  the  differential  path  length  along  the  field  line 
connecting  A  to  C,  B,  is  the  field  strength  along  the  magnetic  field  line,  and  B  is  the  field 
strength  at  A  and  G  Mcllwain  showed  that,  for  a  dipole  field  of  magnetic  moment  M,  the 
required  equatorial  intersection  distance  r  is  given  by  a  relation  of  the  form: 

r*B/M  -  f(I3B/M) 

Although  the  function  f  cannot  not  be  expressed  exactly  in  analytic  form,  excellent 
approximations  are  available  [Mcllwain,  1961,  1966;  Hilton,  1971],  By  Mcllwain’s 
definition,  the  functional  relationship  between  L,  B,  and  I  in  the  given  field  is  the  same  as 
that  between  r,  B,  and  I  in  the  dipole:  Thus,  for  a  specified  magnetic  moment  M,  L  is  a 
function  of  the  point  values  B  and  I. 
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In  the  definition  of  L  presented  here,  we  have  deliberately  omitted  a  specification  of  the 
reference  dipole  moment  M  to  emphasize  the  present  ambiguity  in  this  regard.  Mcllwaift 
initially  selected  a  value  that  was  valid  in  I960.  This  fixed  value  was  imbedded  in  software 
that  he  provided  to  other  users.  Meanwhile  other  software  has  been  distributed  which  uses 
the  updated  dipole  moment,  derived  from  the  dipole  terms  of  the  model  used  to  compute  the 
field.  ‘Hie  consequences  of  ambiguity  on  the  computed  value  of  L  will  be  discussed  in  detail 
in  a  later  sub-section. 


•1.2.1  Elaboration 

Contrary  to  what  we  have  just  said,  let  us  assume  that  we  have  a  universally  accepted 
definition  if  M.  Then  let  us  clarify  the  rest  of  the  definition  of  L  by  considering  two 
magnetic  field  models:  the  model  actually  specifying  tiie  magnetic  field  vector  at  the  point  in 
question,  and  a  reference  dipole  model  with  magnetic  moment  M.  For  either  model,  vve  can 
compute  the  field  intensity  at  any  point,  except  at  the  origin  where  it  becomes  infinite. 

We  can  also  compute  the  invariant  integral  at  any  point  lying  on  a  closed  field  line.  In  the 
dipole  field  all  lines  are  closed  except  the  one  passing  through  the  poles.  On  each  closed 
dipole  line  the  field  strength  minimizes  at  one  point,  on  the  equator,  about  which  the 
original  point  and  its  conjugate  are  symmetrically  located.  In  a  more  general,  non-dipole 
field,  there  may  be  many  open  lines.  For  these  we  cannot  define  either  1  or  L,  so  we  will  not 
concern  ourselves  with  them.  Open  lines  cannot  occur  in  the  radiation  belt,  for  particles 
would  escape  along  them.  Most  field  lines  in  the  radiation  belts  contain  only  one  point  of 
minimum  field  strength,  not  necessarily  on  the  dipole  equator,  in  exceptional  cases,  near  the 
surface  of  the  earth,  there  may  be  more.  In  any  event,  we  may  trace  the  field  line  from  our 
original  point,  in  the  direction  of  decreasing  field  strength,  past  one  or  more  points  of 
minimum  strengm,  and  on  to  the  conjugate  point,  live  first  point  encountered  where  the  field 
strength  matches  that  at  the  original  point.  '12160  we  compute  the  invariant  integral,  a  line 
integral,  in  accordance  with  the  definition,  between  our  original  point  and  its  conjugate. 

Mclhvnin  showed  that  the  locus  of  points  of  common  B  and  I  in  a  dipole  field  has  a  very 
simple  geometric  property:  it  lies  on  a  shell  of  field  lines  intersecting  the  equator  on  a  circle 
of  radius 


r  -  [(M/B)f(I3B/M)J 1/3 

The  value  of  L  of  our  point  of  interest,  in  the  actual  field,  is  set  equal  to  this  radius,  in  RE, 
of  the  shell  occupied  by  the  points  of  the  same  B  and  I  in  the  reference  dipole  field.  Our 
point  does  not  actually  lie  on  this  shell,  unless  the  actual  field  coincides  with  the  reference 
dipole  field.  However,  in  many  instances,  the  radial  distance  of  the  actual  equator  crossing 
of  the  field  line  is  close  to  L 

For  self-conjugate  points,  1-0,  L  is  easily  obtained  from  B.  In  a  dipole  field,  the  self¬ 
conjugate  points  are  on  the  equator,  where 

B  »  M/r3 
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Therefore,  the  equatorial  radius  of  the  dipole  shell  containing  the  1-0  points  of  specified  B 

is 


r  -  (M3J),fl 

We  therefore  must  have  for  1-0  points  in  our  actual  field: 

L  - 

Incidentally  it  folk)**  from  these  considerations  that 
f(0)  -  1 


4.1?.  Physical  Significance 

'Hie  L  parameter  is  related  to  the  properties  of  the  adiabatic  motion  of  trapped  particles  in  u 
sufficiently  slowty  varying  (both  spatially  and  temporally)  magnetic  field.  Under  these 
conditions  the  motion  of  the  particle  can  be  decomposed  into  three  components:  the 
promotion  around  the  magnetic  field  lines,  the  bounce  motion  between  the  conjugate  mirror 
(Xiints,  and  the  longitudinal  drift  motion  around  the  earth.  Associated  with  each  of  these 
components  of  the  motion  is  an  adiabatic  invariant  which  is  approximately  a  constant,  ’the 
first  two  invariants,  associated,  respectively,  with  the  gyromotion  and  the  bounce  motion,  are: 

p  -  p7sinJa/(2mB) 

J  -  2plw 

Here  p  is  the  particle's  momentum,  a  is  the  pitch  angle  (angle  between  the  velocity  and  the 
magnetic  field  vector),  m  is  the  mass,  B  is  the  magnetic  field  strength,  and  Im  is  the  invariant 
integral  at  the  particle’s  mirror  points,  where  the  particle's  pitch  angle  is  90*  and  the 
magnetic  field  strength  Bm  is  given  by  conservation  of  the  first  adiabatic  invariant  p: 

Bw  -  B/sin^. 

In  the  absence  of  an  electric  field,  the  magnetic  field  is  static  and  the  momentum  p  is  a 
constant  of  (he  motion.  Then  the  particle's  mirror  point  magnetic  field  and  invariant 
integral,  Bm  and  Im,  are  also  constants  of  the  motion,  from  which  it  follows  that  its  mirror 
point  L  value,  L*,  is  also  a  constant.  Therefore  the  particle's  motion  is  confined  to  a  "shell" 
of  magnetic  field  lines  which  contain  the  locus  of  points  of  the  specified  Bm  and  L^. 

In  a  pure  dipole  field  of  magnetic  moment  M,  L^,  is  the  equatorial  radius  of  this  shell,  in 
earth  radii,  since  the  field  lines  intersect  the  equator  at  this  distance.  This  simple  geometric 
interpretation  is  the  reason  that  L,,  has  long  been  favored  over  Im  as  a  parameter  for 
modeling  the  radiation  belts.  In  this  dipole  field,  all  particles  found  at  a  point  possess  the 
same  Lm,  since  L  is  constant  along  the  field  line.  In  a  non-dipole  field,  however,  the  L 
parameter  is  not  in  general  equal  to  the  equatorial  intersection  distance  of  the  field  line 
through  the  point,  nor  is  L  a  constant  along  the  field  line,  although  both  of  these  properties 
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,irc  often  approximately  satisfied  Therefore  tkc  L  parameter  at  a  point  u  the  conferred 
[,a  only  for  tkc  particle*  mirroring  at  that  point,  that  it,  tkotc  parlictct  with  90*  pitch 
angle*  at  that  point.  Tkc  appropriate  lm  parameter  value  for  particle t  with  tmaller 
pitch  angle*  muit  he  computed  at  their  particular  mirror  point*,  and  in  general  mill 
therefore  he  different  from  the  L  parameter  value  at  the  point  of  ohicrvation. 


4.3  computational  methods 

From  the  definition  of  L  given  here,  it  is  evident  that  four  steps  are  required,  in  the 
M lowing  order: 

1.  Compute  the  model  field  vector  and  magnitude  B  at  the  point  in  question. 

2.  Trace  the  field  line  through  the  point  to  its  conjugate  point,  storing  the  magnetic  field 
magnitude  B,  at  points  along  this  path. 

3.  Compute  the  line  integral  along  the  field  line  to  obtain  I. 

4.  Compute  L  from  I,  B,  and  M. 

Step  1  is  reasonably  straightforward,  given  the  mathematical  specification  of  the  field  model, 
and  will  not  be  discussed  further  here.  Step  2,  the  field  line  tracing,  is  the  most  time 
consuming  computationally,  since  it  requires  solution  of  a  vector  system  of  nonlinear 
ordinary  differential  equations  expressing  tire  local  tangency  of  the  magnetic  field  vector. 
Step  3  requires  only  a  quadrature  integration  along  the  field  line,  having  determined  from 
step  2  the  field  strength  as  a  function  of  the  distance  along  the  field  line  from  starting  point. 
Finally,  in  step  4,  we  must  compute  L,  given  I,  B,  and  M.  The  shortest  of  the  four  s»eps  of 
our  procedure,  this  final  step  requires  selection  of  a  algorithm  for  computing  I.  with 
sufficient  accuracy,  and  tpccijiculion  of  the  magnetic  moment  M  of  the  reference 
dipole  field. 

The  expense  of  the  computation  of  L  directly  from  its  definition  may  be  avoided  by 
interpolation  of  prestored  tables,  if  the  same  model  will  be  employed  many  times  more  than 
the  number  of  tabular  entries  required  ro  interpolate  L  to  sufficient  accuracy. 


4.3.1  Step  2  -  Field  One  Tracing 

'Hie  desired  field  line  is  tire  solution  of  the  ordinary  differential  equations 
dx/ds  -  B/B 
dy/ds  -  B/B 
dz/ds  “  Bj/B 

satisfying  the  initial  conditions  that  x,  y,  z  at  s-0  are  the  coordinates  of  the  point  at  whicn 
we  wish  to  compute  L  Here  x,  y,  and  z  are  the  cartesian  cooidinates  of  a  point  on  the  field 
line,  and  s  is  the  distance  of  the  point,  measured  along  the  field  line,  from  our  original  given 
point,  with  a  positive  sign  if  the  tracing  is  parallel  to  the  magnetic  field  vector,  a  negative 
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sign  if  the  tracing  is  antiparallcl.  Thus,  the  field  line  is  defined  as  a  series  of  points,  specified 
by  their  Cartesian  coordinates,  as  functions  of  the  independent  variable  s. 

4.3. 1.1  Kluge's  method 

Programs  OPTRACE  and  HMIN  solve  these  equations  as  indicated.  The  principal 
difference  between  the  two  programs  arises  because  OPTRACE  must  obtain  the  100  km 
intercepts  of  the  field  line  as  part  of  the  CRRES  requirements.  Therefore  it  must  trace  the 
field  all  the  way  to  these  intercepts,  in  both  hemispheres.  HMIN,  concerned  only  with 
computing  L,  traces  only  to  the  conjugate  point. 

SI  1ELLC,  using  Kluge’s  method,  invokes  a  change  of  coordinates  to  the  inverse  variables. 


£  "  x/r2 
■n  -  y/r2 
{  -  z/r2 

Here  x,  y,  and  z  arc  in  a  dipole-oriented  system,  the  z  axis  parallel  to  the  dipole.  From  the 
inverse  coordinates  we  cm  readily  recover  the  direct  coordinates: 


x  "  C/p2 


y  -  ttf)2 
z  -  {/p2 


where 


p2  -  f}  +  -rj2  +  £2  «  1/r2 

‘lliere  is  also  a  change  in  the  independent  variable  from  s  to  the  variable  t,  such  that 
dt  -  p5dt'  ■  psd$/B 

’Ihe  differential  equation  for  $  becomes 
d£/dt  -  3W/a$  +  $W/p2 


where 


W  -  -rV 


and  V  is  the  scalar  potential  from  which  the  field  is  derived: 


B  «  -VV 


Fully  analogous  equations  follow  for  q  and  5.  If  the  field  is  not  derivable  from  a  scalar 
potential,  as  when  the  external  pc  ;ion  (due  to  ring  currents,  etc.)  is  included,  we  instead 
would  obtain 
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d$/dt  -  p-3((ae/ax)B,  +  (ae/ay)By  +  (a^az)Bj 

At  (his  point  the  final  change  of  coordinates  is  made,  replacing  £  and  t]  by 


u  - 


v  ■  iyp,/2 


In  n  pure  dipole  field,  u  and  v  are  constants  along  the  field  line: 


u  -  cos<KL1/J 


v  m  s\n$flJn 


Also,  d£/dt  is  shown  to  be  positive  everywhere  (at  least  for  the  pure  dipole  case).  Since  (  and 
t  are  therefore  single-valued  functions  of  each  other,  (  can  serve  as  the  independent  variable 
instead  of  t.  Then  from  the  definitions  of  u  and  v,  equations  are  obtained  for  du/d(  and 
dv/d{: 


du/d(  -  [p-,/2d£/dt  -  (l/2)p^dpWt]/(d{«t) 
dv/d$  -  [p->«dVdt  -  (l^-^Tid^dtKd^dt) 

The  right-hand  sides  of  these  equations  can  be  expressed  completely  in  terms  of  u,  v,  and  (. 
First  note  that 

dp/dt  ■*  ((d(/dt  +  ndr/dt  +  (d(/dt)/p 

'Hie  derivatives  d£/dt,  dVdt,  and  d£/dt  are  given  above  in  terms  of  p,  (,-n,  and  (.  Thus  du/d£ 
and  dv/d(  are  expressible  entirely  in  terms  of  these  variables.  To  eliminate  p,  £,  and  -q,  make 
the  substitutions: 

£  -  up>/2 

m  yp,/J 

P  -  (l/2)(u2  +  v2  +  [(u2  +  v2)2  +  4£2]1/2) 

The  expressions  for  £  and  n  come  from  inversion  of  the  definitions  for  u  and  y.  The 
expression  for  p  comes  from  solving  the  quadratic  equation 

p2  -  {2  +  H2  +  (2  -  p(u2  +  v2)  +  (2 

With  these  substitutions  we  obtain  ordinary  differential  equations  for  u  and  v  in  terms  of 
themselves  and  the  independent  variable  (.  Solution  of  these  gives  u(()  and  v(().  As  stated 
before,  u  and  v  are  constant  along  a  dipole  field  line.  In  the  real  world,  we  expect  them  to 
be  sufficiently  slowly  varying  that  very  large  steps  can  be  used  in  the  tracing.  Typically  10 
steps  seem  to  be  sufficient  for  this  method,  while  for  the  direct  method  used  in  OPTRACE 
and  HMIN,  the  number  of  steps  is  ~100.  Note  in  addition  that  only  two  equations  need  be 
solved  here,  while  in  the  direct  approach  we  are  solving  three. 
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•0.1.2  Method!  of  numerical  integration 
'Hie  equations  to  be  solved  are  of  the  form 
dy/dx  -  g(y,x) 

The  principal  feature  of  these  equations  is  that  the  right-hand  side,  the  derivative  of  the 
solution,  depends  on  the  value  of  the  solution  itself.  A  second  feature  relevant  to  our 
particular  problem  is  that  the  derivatives  (right-hand  sides)  require  the  time-consuming, 
evaluation  of  the  magnetic  field. 


OPTRACE/HMIN 

OPTRACE  and  HMIN  i jt  the  open  Adams  4th  order  numerical  integration  method  [Prett, 
el.  ai,  1986;  Hildebrand,  1956].  This  method  generates  a  solution  on  a  sequence  of 
uniformly  spaced  points.  In  each  step  the  solution  at  a  new  point  is  generated  from  the 
value  of  the  solution  at  .he  previous  point  and  the  first  derivatives  of  the  solution  at  the 
previous  four  points.  Only  one  derivative  evaluation  is  required  per  step. 

'Hie  method  is  called  "open"  to  distinguish  it  from  other  methods  which  are  "dosed",  since 
they  require  also  the  derivative  of  the  solution  at  the  new  point.  However  the  first  derivative 
of  the  solution  value  at  any  point  can  be  obtained  only  if  the  solution  value  itself  is  known, 
since  (he  right-hand  sides  of  the  equations  are  functions  of  the  solution  value  as  well  as  the 
independent  variable.  The  dosed  methods  are  often  used  in  what  are  called  predictor- 
corrector  schemes,  in  which  an  open  method  is  used  to  "predict"  the  solution  at  the  new 
point.  Then  the  predicted  solution  value  at  the  new  point  is  used  to  estimate  the  derivative 
there  for  use  in  the  dosed  formula. 

The  method  is  called  4th  order  to  indicate  that  it  is  accurate  through  the  4th  power  in  the 
step  size,  i.e.,  its  error  is  proportional  to  the  5th  power  in  the  step  size. 

Since  the  solution  must  be  known  at  4  previous  points  in  order  to  obtain  it  at  a  new  point,  it 
must  be  generated  at  the  first  4  points  in  the  sequence  by  another  method.  Actually,  the 
solution  is  given  to  us  at.  the  first  point,  the  boundary  value.  Thus  we  must  generate  from  it 
the  solution  values  at  the  next  three  points.  In  OPTRACE  and  HMIN  this  is  done  by  the 
classical  Runge-Kutta  4th  order  method  [Prett,  el.  at.,  1986;  Hildebrand,  1956],  which 
allows  one  to  obtain  the  solution  ai  a  new  point  from  just  its  value  and  derivative  at  the 
previous  point.  'Hie  process  requires  evaluations  of  the  derivatives  for  two  estimated  solution 
values  at  the  midpoint  of  the  step,  plus  one  evaluation  for  an  estimated  solution  value  at  the 
end  of  the  step,  in  addition  to  evaluation  with  the  known  solution  value  at  the  beginning 
point  of  the  step.  Thus  each  step  requires  four  derivative  evaluations,  compared  to  just  one 
per  step  for  the  open  Adams  method.  This  makes  Runge-Kutta  considerably  slower  than 
Adams,  but  it  has  the  flexibilty  of  not  requiring  a  previous  multiple-point  history. 

The  step  size  h  is  given  by: 

h  -  FLU[SEGS  x  f(FLL)] 


where 
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FLL  «  estimated  field  line  length  (Re)  between  100  kin  intercepts 
SEGS  "  number  of  steps  desired  when  FLL  *  2  RE 
f(FLL)  -  (-273) FLL  +  7/3,  FLL  <;2  RE 
-  (1/18)FLL  +  BJ9,  FLL>2RE 

‘Hie  expression  for  h  reflects  the  finding  [Radex,  Inc.,  1987]  that  uniform  accuracy  was 
obtained  with  a  number  of  steps  that  minimizes  for  FLL  ■  2  and  increases  linearly  to 
approximately  double  that  amount  at  the  shortest  and  longest  lengths  to  be  considered  for 
CRRES.  'litis  shape  function  is  reflected  by  the  function  f(FLL),  while  SEGS,  nominally  set 
to  100,  but  easily  adjusted  in  the  code,  gists  the  desired  number  of  steps  at  FLL  -  2  RE. 

'ihe  procedure  for  estimating  the  field  length  depends  on  the  dipole  L  parameter,  Lj,  of  the 
initial  point  [McNeil,  1986].' 

Lj  -  r/sin^fl 

where  r  is  the  radial  distance,  and  0  is  the  magnetic  colatitude.  For  Lj  £  2  RE,  FLL  was 
found  empirically  to  fit  veiy  well  the  estimate 

FLL  -  2.77  Lj  -1.86 

For  lesser  values  the  situation  is  quite  complicated,  due  to  the  varied  effects  of  the  non-dipole 
terms  in  the  field.  In  this  case  FLL  was  fit  to  parabolic  functions  of  Lj,  with  coefficients 
tabulated  on  a  7  x  13  magnetic  latitude-longitude  grid.  FLL  at  any  latitude  and  longitude  is 
computed  by  bi-lincar  interpolation  of  the  four  nearest  neighbors  on  the  grid. 

SHELLC 

The  third  order  Adams  predictor-corrector  method  is  used  to  obtain  the  solution  at  a  new 
point,  given  the  solution  at  the  previous  point  and  the  derivatives  at  the  previous  three.  'Hie 
process  is  started  as  follows: 

1)  the  solution  is  estimated  1/2  step  forward  of  the  starting  point  by  linear  extrapolation, 
using  the  value  and  derivative  at  the  starting  point; 

2)  The  solution  is  then  obtained  one  step  backward  from  the  starting  point,  by  quadratic 
extrapolation  using  the  starting  point  solution  value  and  derivative,  and  the  estimated 
derivative  at  the  point  1/2  step  forward; 

3)  The  solution  one  step  forward  of  the  starting  point  is  obtained  by  cubic  extrapolation 
using  the  solution  value  at  the  starting  point  and  the  derivatives  at  all  three  previously  known 
points.  The  point  1/2  step  forward  is  now  discarded,  and  the  remaining  three  arc  available 
to  begin  the  predictor-corrector  process. 

'Hie  step  size  for  all  runs  was  0.2  Rg1,  the  value  set  in  the  code  we  received.  Since  this 
probably  was  arrived  at  after  extensive  testing,  we  saw  no  reason  to  change  it. 
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<1.3.2  Step  3  -  Compulation  of  I 


OPTRACE/HM1N 

’llic  interval  between  the  starting  point  and  its  conjugate  point  is  divided  into  two  regions: 

1.  the  largest  subinterval  beginning  at  the  starting  point  which  contains  an  even  number  of 
steps; 

2.  the  subinterval  beginning  at  the  end  of  region  1  and  ending  at  the  conjugate  point. 

'Hie  conjugate  seldom  coincides  with  any  grid  point;  it  is  found  by  quadratic  fit  of  the  field 
strength  B,  through  the  three  closest  points,  and  solving  for  s  such  that  B,  «  B. 

The  integral  through  region  1  is  performed  by  Simpson’s  Rule  [Prat,  tt.  ai,  1986; 
Hildebrand,  1956],  which  is  a  two-step  formula.  The  midpoint  of  region  2  is  located  by  a 
single  application  of  Runge-Kutta;  then  one  application  of  Simpson’s  Rule  computes  the 
integral  through  this  portion. 

If  the  conjugate  point  lies  within  the  first  two  steps  a  different  procedure  is  followed.  'lire 
quantity  l-B/B  is  fit  to  a  parabola,  from  which  the  location  and  value  of  the  minimum 
(Bm[n)  are  determined.  'Hie  invariant  integral  is  then  given  analytically  by 

I  -  (ir/4)t(l-Bm|n/B),/2 

where  t  is  the  ardength  from  the  starting  point  to  the  conjugate. 

S  HELLO 

In  terms  of  the  transformed  variables  the  invariant  integral  is  written: 

i- j(i-Bt/B)w  mwmm 

where 

Q  -  [(d£/dt)2  +  (drydt)2  +  (dC/dt)2p 

It  is  shown  that  the  field  strength  B(  at  a  point  on  the  field  line  am  be  written  as 
B,-p?Q 

The  step  size  used  for  the  tracing  is  too  large  for  the  computation  of  I,  because  of  the 
rapidly  varying  quantities  in  the  integrand.  To  avoid  having  to  trace  at  a  o. nailer  step  size, 
which  would  defeat  the  purpose  of  the  algorithm,  Kluge  defines  three  slowly-varying 
interpolation  functions: 

C  *  u2  +  v2 

D  «  Q/[l  +  3((/p)2p 

E  “  D(1  +  £2/p2)/(d£/dt) 
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In  a  dipole  field,  C  becomes  1/L,  imd  D  and  E  both  become  equal  lo  the  dipole  moment, 
llte  inverse  radius  p  has  previously  been  given  in  terms  of  the  slowly-vatying  u,  v,  and  the 
independent  variable  {.  We  an  re-express  this  result  here  as: 

p  -  (1/2)(C  +  [C2  +  4{2]>'2) 

Huts  we  an  compute  p  analytially  at  any  pointC  once  we  know  C  Next,  Q  an  be  obtained 
analytically  from  D,  p,  and  (;  d{/dt  an  be  computed  from  D,  E,p,  and  and  13,  an  be 
computed  from  p  and  Q.  Thus  the  entire  integrand  an  be  computed  analytically  at  any 
point  {  from  these  interpolation  functions.  Four-point  interpolation  is  used  for  C,  three- 
point  interpolation  for  D  and  E 

The  integration  is  performed  by  the  trapezoidal  rule  [Pro*,  «t.  a/.,  19S6).  Since  this 
method  performs  the  integration  one  step  at  a  time,  the  step  size  an  be  varied  conveniently. 
It  is  chosen  to  be  kp,  where  k  Is  presently  set  to  0.03. 

"Hie  beginning  and  ending  intervals  are  handled  by  a  slight  modifiation  to  account  for  the 
square  root  singularity  in  integrand  -  it  approaches  zero  as  the  square  root  of  the  distance 
from  the  end  point.  If  xt  is  the  length  of  the  interval,  and  x  is  the  distance  from  the 
endpoint  where  the  integrand  vanishes,  then  the  integrand  is  approximated  by 

G  -  a[x/x,]w 

where  a  is  the  value  of  the  integrand  at  the  other  end  of  the  interval.  ’Ilien  the  integral  over 
this  portion  is 


i  -  (2/3)ax, 

If  the  conjugate  point  is  within  one  integration  step  of  the  starting  point,  I  is  set  to  zero. 


•1.3.3  Step  4  -  Determination  of  L  from  B  and  I 

Although  an  exact  matliematial  expression  for  L  is  not  available,  excellent  approximations 
have  been  provided  by  Mcllwain  [1961,  i966]  and  Hilton  (1971).  All  three  codes  employ 
Mcllwain’s  1966  version: 

Y  -  ln(L3B/M-l)  -  2n.0’  anX" 


where 


X  -  ln(I3B/M) 

and  the  coefficients  a„  are  given  in  the  table  reproduced  here  from  Mcllwain  [1966].  'lire 
last  row  in  the  table  indiates  the  maximum  error  in  Y.  These  numbers  translate  into  a 
maximum  relative  error  (AL/L)  of  ~1.6xl(H. 

Originally,  the  somewhat  simpler  formula  given  by  Hilton  [1971], 

L3B/M  -  1  +  a,Xw  +  ajX2*  +  a3X 


where 


35 


X  -  PB/M 
a,  -  1.35CM7 
a2  -  0.465376 
a3  -  0.0475455 

was  employed  In  OPTRACE  and  HMIN.  Table  1  of  Hilton's  paper  quotes  maximum  error 
AIVL  — 10*4.  Hotter,  comparisons  with  results  obtained  using  Mcllwain's  formula  revealed 
relative  differences  ~  2xl0*3.  'ITtus  we  reverted  to  the  tried  and  true  Mcllwain  formula. 
Recently  we  have  discovered  that  the  differences  may  due  to  an  error  in  a2  in  the  computer 
code  of  Hilton's  formula.  Although  further  testing  may  warrant  reinstatement  of  the  Hilton 
formula,  the  saving  in  computation  is  minimal,  since  this  step  requires  a  very  small  fraction 
of  the  total  effort  required  to  compute  L  Therefore  we  have  for  the  time  being  retained 
Mcllwain's  formula. 

4.3.3.1  The  disappearing  dipole 

Historically  live  constancy  of  the  Earth’s  magnetic  dipole  has  been  taken  for  granted.  As  in 
human  relations,  this  may  prove  to  be  fatal,  for  the  dipole  moment 

M  -  ((g,0)2  +  (g,l)2  +  (h,')2)w 

has  been  decreasing  at  the  rate  of  approximately  0.05  %  per  year.  The  predicted  1990  value, 
computed  from  the  dipole  terms  of  the  extrapolated  IGRF  1985  model,  is  30299  (nT-Rc3), 
while  the  value  used  by  Mcllwain,  in  these  units,  is  31165.3.  This  latter  value  was  imbedded 
into  Mcllwain's  subroutine  CARMEL,  which  be  has  over  the  years  provided  to  many  users. 
Thus  this  constant  has  become  a  fixture  in  many  L  computation  codes,  including  the 
SHELLC  package  sent  to  us.  In  OPTRACE  and  HMIN,  on  the  other  hand,  we  have  used 
the  updated  value  for  M,  given  above,  as  recommended  by  Hilton  (1971).  Since  we  have  only 
recently  discovered  that  it  was  Mcllwain's  old  M  value  that  was  used  in  SHELLC  many  of 
the  results  of  that  code  presented  in  this  report  were  computed  with  this  value. 
Consequently  part  of  the  differences  between  SHELLC  and  OPTRACE/HMIN  results  are 
undoubtedly  due  to  this  difference  in  the  M  values  used.  We  have  since  modified  SHELLC 
to  use  tire  updated  dipole  moment  as  in  OPTRACE,  and  will  present  comparative  results. 
The  INTEL  interpolation  tables,  undoubtedly,  have  been  generated  using  the  old  dipole 
moment.  Since  we  do  not  presently  have  the  capability  to  regenerate  these,  the  only  INTEL 
results  presented  here  are  with  the  old  dipole  moment. 

It  is  not  totally  obvious  to  us  how  to  select  M.  We  have  given  here  the  two  choices  that  we 
know  have  been  commonly  used.  Mcllwain  [1966]  suggested  that  M  be  chosen  to  minimize 
the  variation  of  L  along  the  field  lines.  Although  attainment  of  a  practical  algorithm  for 
accomplishing  this  in  the  actual  geomagnetic  field  seems  unlikely,  Hilton’s  definition,  used  in 
OPTRACE,  accomplishes  this  in  a  time-varying  pure  dipole  field:  the  L  parameter  at  a  point 
remains  equal  to  the  radial  distance  of  the  equatorial  intersection  of  the  field  line  through 
that  point  as  the  field  changes. 

However,  what  happens  if  we  use  Mcllwain’s  fixed  value,  which  we  may  call  Mm,  instead  of 
the  updated  value,  Mu?  We  have  found  previously,  that  for  NO, 
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L  -  (M/B)«a 


Thus  we  find  that  by  using  Mcllwain's  fixed  value  of  the  dipole  moment,  instead  of  the 
updated  value,  we  obtain  for  points  at  the  equator  (1-0)  L  parameter  values  which  differ 
from  the  radial  distance  by  the  factor  (Mm/M.,)10.  Using  the  values  of  Mu  (1990)  and  Mm 
given  above,  we  find  that  this  factor  is  1  Thus,  in  1990,  use  of  the  fixed  value  Mm  would 
result  in  L  values  nearly  ICS  higher  at  tho  equator  than  those  obtained  using  the  updated 
value  Mu. 

Now  consider  high  latitudes,  in  the  limit  of  large  I  and  B.  Here  both  Mcllwain’s  and 
Hilton's  expansions  lead  us  to 

L3B/M  =  0.0475455 13B/M 
l/I  a  0362273 

'Iltus  L  becomes  independent  of  M  in  the  high  latitude  limit,  while  it  depends  on  M  at  the 
equator.  Therefore  wc  conclude  that  L  is  not  a  constant  along  a  dipole  field  line  unless  the 
value  of  the  dipole  moment  M  employed  in  computation  of  L  from  B  and  I  Is  the  same  as 
the  value  employed  in  computing  B  and  I.  Thus  to  satisfy  Mcllwain’s  criterion  of 
minimizing  the  variation  of  L  along  a  field  line,  it  is  obvious  that  for  a  pure  dipole  field  we 
should  use  the  same  (updated)  value  in  the  final  step  of  obtaining  L  from  B  and  I  as  was 
used  in  computing  the  field.  The  effect  of  using  Mcllwain’s  subroutine  CARMEL,  with  the 
lived  value  Mm  embedded,  in  conjunction  with  model  routines  that  use  the  updated  values  of 
all  coefficients  to  compute  the  field,  is  that  this  condition  would  not  be  satisfied. 

In  a  non-dipole  field,  however,  L  is  not  in  general  constant  along  a  field  line,  even  if 
consistent  values  of  M  are  employed.  For  the  earth’s  internal  field,  however,  L  is 
approximately  constant,  Mcliwain  [1961]  finding  relative  variations  within  1%  along  field 
lines  within  3  Ru  (Hopefully  his  M  value  is  consistent  with  the  1960  model  he  employed). 
However  at  larger  distances,  where  the  field  due  to  external  sources  must  be  included,  the 
variation  of  L  along  a  field  line  an  be  quite  substantial.  For  example,  using  the  Olson- 
Pfitzer  model  for  the  external  field  and  the  IGRF85  for  the  internal  field,  it  was  found,  for  a 
field  line  at  2200  hours  MLT  with  an  840  km  footprint  ut  65°  magnetic  latitude,  that  L  was 
7.64  at  the  footprint  and  8.88  at  the  equatorial  crossing,  for  a  relative  variation  of  more  than 
XlKc.  Therefore  it  isn't  clear  to  us  what  choice  of  M  would  minimize  variations  of  L  along 
field  lines  of  a  realistic  geomagnetic  field  model. 

In  summary,  the  earth’s  dipole  moment  M  has  changed  by  nearly  3%  since  the  time  of 
Mcllwain’s  original  definition.  If  we  use  Mcllwain’s  original  value  in  computing  L  from  B 
and  I,  we  get  near  the  equator  L  values  approximately  1%  higher  than  those  obtained  when 
using  the  updated  dipole  moment.  There  are  several  L  computation  codes  that  use  the  fixed 
Mcliwain  value,  and  some  cocfcs  that  use  the  updated  value.  We  should  therefore  be  warned 
that  the  results  obtained  with  these  codes  may  disagree  due  to  this  disparity.  It  would  be 
wise  to  agree  in  the  future  on  a  standard  choice  to  be  used  by  everyone. 


4.3.4  Interpolation  of  L  -  Program  INTEL 
The  L  parameter  is  expressed  in  the  form 
L  -  d-i[l+f(d,?,4>)] 
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where 


d  -  aJ/(ff2+  ;2)I/2 

«  -  KW*)2  +  (•n-'nP)2J,/2 
cos4>  - 
sin4>  " 

Here  yO  and  ^(0  define  the  open  field  line,  which  goes  over  to  the  dipole  limit  $p  -  Tip 
■  0  at  C-0.  They  are  expressed  as  4th  degree  polynomials  (constant  term  ■  0).  This 
reformulation  removes  the  polar  singularity  from  the  function  which  must  be  interpolated  (i. 
e.,  the  function  f  remains  finite).  The  evaluation  of  f  is  by  quadratic  interpolation  of  a 
fourier  series  in  $  on  a  uniform  £-d  grid.  The  spacing  of  the  grid  is  Ad  ■  ■  0.2,  with 
dwta  -  0.1,  dWM  ■  0.9,  Cmt«  "  *0.9,  -  0.9.  Not  all  points  on  this  grid  are  actually 

employed,  since  some  are  below  the  surface.  A  total  of  700  Fourier  coefficients  are 
tabulated,  10*30  per  point,  at  the  grid  points  above  the  surface. 


38 


5.  L-SHELL  COMPARISONS  OF  THE  MODELS 


The  four  models  (Barraclough  1975,  MAGSAT  1980,  IGRF  1985,  and  Gun  1990)  were 
compared  on  the  basis  of  their  L-shell  evaluation  as  well  as  their  B  field  determination  (see 
Oiapter  3).  Each  of  these  models  was  incorporated  into  software  currently  in  use  at  GL 
To  begin  with,  contour  plots  of  the  differences  between  the  models  were  made  over  a  large 
altitude  range.  Then  line  plots  of  L-shcll  at  selected  longitudes  were  made  for  a  wide  latitude 
range,  ±75°.  To  more  closely  examine  differences  between  the  models  at  latitudes  pertinent 
to  CRRES,  detailed  line  plots  were  made  at  five  altitudes  for  several  longitudes  with  latitude 
only  spanning  ±20°.  In  the  difference  comparisons  mentioned  previously,  IGRF  19.85  was 
used  as  the  standard,  'nuts,  contour  plots  of  the  IGRF  19.85  L*siiell  over  the  whole  latitude 
and  longitude  range  were  done.  These  were  made  for  1980,  1985,  and  1990  to  look  for 
temporal  variations  in  L 


5.1  OVERALL  MODEL  DIFFERENCES 

Figures  5.1a  and  5.1b  show  the  differences  between  the  models  for  1990  at  350  km  and 
1U.U00  km,  respectively.  IGRF  19S5  was  used  as  the  reference  because  it  is  the  most  widely 
used  model  and  was  adopted  as  the  standard  by  the  IAGA  Working  Group  1.  'Hie  contour 
plot  in  the  lower  right  hand  panel  shows  L-shell  as  determined  by  IGRF  1985.  The 
remaining  three  panels  show  the  differences  between  the  other  three  models  and  this  IGRF. 

'Hie  upper  left  panel  shows  the  absolute  differences  between  Darraclough  and  IGRF.  'Hi is 
model  differs  the  most.  Even  so,  these  differences  tend  to  be  less  than  a  percent  or  two 
from  the  standard.  The  agreement  deteriorates  as  one  approaches  the  poles.  This  is  what  is 
expected  since  L  increases  toward  the  poles  and  the  models  themselves  are  not  as  good 
nearer  the  poles  (sec  Chapter  2).  Whether  Barraclough  is  higher  or  lower  than  IGRF 
depends  upon  the  hemisphere.  This  may  be  due  in  part  to  the  change  in  the  position  of  the 
offset  dipole.  'Hie  dipole  moment  and  the  displacement  of  the  offset  dipole  are  known  to  be 
changing  with  time  (sec  sections  4.3.3.1)  (Pinto,  el.  al.,  1989).  For  this  comparison,  the 
software  determines  the  position  of  the  offset  dipole  and  the  dipole  moment  from  the 
coefficients.  Thus,  different  offsets  would  be  expected  to  lead  to  systematic  differences  in  L 
by  hemisphere.  Merc,  Barraclough  is  too  low  in  the  North-East  and  South-West  quadrants. 
It  is  too  high  in  the  North-West  and  South-East  quadrants.  Note,  these  differences  are 
prima.ily  attributable  to  differences  between  the  models  themselves.  Tlie  differences  in  the 
offsets  add  to  these  larger  differences. 

The  upper  right  panel  shows  differences  between  MAGSAT  and  IGRF.  Here  the  agreement 
is  seen  to  be  a  little  better  than  that  with  Barraclough.  More  area  is  covered  by  the  +.001 
contours  and  the  larger  contours  are  at  higher  latitudes.  Still,  the  percent  difference  from 
the  standard  is  about  the  same:  generally  less  than  one  percent.  Again  the  sign  of  the 
difference  depends  on  the  hemisphere  with  MAGSAT  too  low  in  the  North-E^st  and  South- 
West  quadrants  and  too  high  in  the  North-West  and  South-East  quadrants  (roughly). 

Finally,  the  best  agreement  is  seen  between  Giin  and  IGRF  in  the  lower  left  panel  of  Figure 
5.1a.  For  the  most  part  the  absolute  differences  are  less  than  .001  with  higher  differences 
very  close  to  the  poles.  This  is  an  order  of  magnitude  better  than  the  two  previous  models 
discussed.  The  relationship  between  the  sign  of  the  difference  and  geographic  position  is  not 
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LONGITUDE  (DEGREES)  LONGITUDE  (DEGREES) 

Figure  5.1a.  Contour  plots  of  model  L-sltcll  differences  in  Earth  radii  (R^)  covering  the  full  latitude  and 
longitude  ranges.  IGRF  19S5  is  used  as  a  reference  model.  The  comparisons  arc  done  for  an  altitude  of 
350  km  (perigee)  at  1990.0. 
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LONGITUDE  (DEGREES)  LONGITUOE  (0  EG  REES) 

Figure  5.1b.  Contour  plots  of  model  L-shcll  differences  in  Earth  radii  (R^  covering  the  full  latitude  and 
longitude  ranges.  IGRF  19S5  is  used  as  a  reference  model.  The  comparisons  arc  done  for  an  altitude  of 


so  clear  cut  as  previously.  These  differences  are  generally  less  than  a  tenth  of  a  percent  of 
the  IGRF  values.  The  instruments  on  CRRES  will  not  be  able  to  discem  L-shell  differences 
at  this  order  of  magnitude. 

Figure  5.1b  is  similar  to  5.1a,  tbo  difference  is  in  the  altitude  at  which  the  comparisons  are 
made.  Mere,  one  secs  the  comparisons  done  for  10,000  km.  Again,  IGRF  1985  is  the 
standard  with  the  other  throe  ntodefe  compared  to  it.  The  L-shcll  contours  obtained  from 
IGRF  1965  are  shown  in  the  lower  right  panel.  The  absolute  differences  contoured  in  the 
other  three  plots  are  all  smoother  than  those  seen  at  350  km.  At  higher  altitudes  the 
contributions  to  the  internal  field  from  the  crustal  components  are  negligible  in  comparison 
to  the  dipole  component,  thus  the  smaller  scale  detail  washes  out  and  the  contours  become 
smoother.  Tire  basic  distribution  of  the  contours  has  remained  the  same  between  350  km 
and  10,000  km  (i.e.  the  distribution  of  tire  positive  and  negative  differences  is  still  correlated 
with  the  offset  dipole).  However,  in  all  three  plots,  the  differences  ,jve  increased.  L 
increases  with  altitude,  so  even  though  the  absolute  differences  have  increased  somewhat,  the 
percent  differences  arc  about  the  same.  Thus,  one  again  sees  that  for  CRRES,  differences 
between  Cain  and  IGRF  will  not  be  discernible. 


5.2  COMPARISONS  AT  SPECIFIED  LONGITUDES 

Four  line  plots  at  selected  longitudes  were  made  of  L-sliell  versus  latitude  (see  Figure  5.2). 
At  this  sale,  the  agreement  between  the  models  is  so  close,  that  only  one  line  is  apparent 
even  though  all  four  models  have  been  plotted  on  each  plot.  These  plots  were  made  using 
the  high  latitude  data,  thus,  one  an  see  longitudinal  variations  in  L  ‘Hie  magnetic  poles  are 
near  the  two  longitude  selected  for  the  right  panels.  Hence,  a  sharp  increase  in  L  is  seen  in 
the  southern  hemisphere  in  the  upper  plot  and  in  the  northern  hemisphere  in  the  lower  plot. 
The  two  left  panels  are  not  near  the  magnetic  poles  in  longitude  so  only  L-shells  £  10  are 
seen. 

Line  plots  were  also  made  at  a  smaller  scale  to  look  at  differences  between  ±20°  (see  Figure 
5.3a  -  5.3c).  Figure  5.3a  shows  L-shell  versus  latitude  at  350  km  and  80°  east  longitude  for 
1990  and  19S0  (upper  and  lower  panels,  respectively)  for  all  four  models.  In  i960  (lower 
panel),  the  agreement  at  (his  sale  is  so  good  that  only  one  line  is  apparent.  This  agreement 
deteriorates  somewhat  by  1990  (upper  panel)  as  is  expected  as  one  gets  farther  away  from 
the  epochs  for  which  the  models  are  derived  and  differences  in  their  secular  variation  are 
manifested.  Nonetheless,  differences  can  only  be  seen  between  -20°  and  4°.  Here, 
Barradough  has  the  highest  L  values  (solid  line),  MAGSAT  is  slightly  lower  (dotted  line), 
and  Cain  and  IGRF  are  lower  still  (dashed  and  dash-dot  lines  plotted  on  top  of  each  other), 
'litis  is  what  is  expected  with  a  decreasing  dipole  moment.  A  higher  moment  will  give 
higher  L  values.  Since,  Barradough  is  the  oldest  model  it  is  expected  to  give  the  highest  L 
values.  Cain  and  IGRF  are  based  on  very  similar  data  sets  (see  Chapter  2).  Even  though 
Otin’s  model  incorporates  a  discontinuity  in  the  secular  variation  after  1983,  the  shift  in  L 
due  to  this  continuity  is  not  yet  apparent  by  1990. 

Figure  5.3b  is  similar  to  Figure  5.3a,  the  difference  is  in  altitude.  'Fills  shows  the  models  at 
10,000  km.  Again,  there  is  no  discernible  difference  between  the  models  at  19S0.  They  all 
plot  on  top  of  each  other.  By  1990,  differences  again  appear  between  -20°  and  4°  latitude. 
Barradough  gives  the  highest  L  values  with  MAGSAT  lower  and  Cain  and  IGRF  lower  still. 
L  increases  with  altitude,  thus  at  10,000  km,  L  is  nearly  2.0  Rg  higher  than  at  350  km. 
However,  the  differences  between  the  models  are  about  the  same  as  those  at  the  lower 
altitude.  Cain  and  IGRF  still  overplot  each  other,  so  even  at  higher  altitudes  one  annot  see 
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Figure  5.3a.  Line  plots  similar  to  those  in  Figure  5.2,  but  over  x  smaller  latitude  range.  Even  at  this 
rcsolution^the  differences  between  the  models  are  very  small.  The  top  panel  shows  the  differences  at  350 
km  and  80*  east  longitude  for  1990.0.  The  bottom  panel  is  the  same  plot  done  for  epoch  19S0.0. 
Discrepancies  arise  over  time,  but  the  differences,  where  distinguishable,  are  still  quite  small. 
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Figure  5.3b.  Line  plots  very  similar  to  those  in  Figure  5.3a,  but  for  an  altitude  of  10.000  km  rather  than 
iiO  km.  The  top  panel  show  the  differences  seen  at  1990  and  the  bottom  panel  show  differences  at  19S0 
for  S0°  east  longitude.  Again,  differences  increase  with  time,  but  not  by  much.  The  altitude  docs  not  seem 
to  affect  the  differences  much;  the  same  relationship  between  the  two  epochs  is  seen  for  10,000  km  as  for 
350  km. 
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the  effects  of  the  jerk  incorporated  in  Cain’s  model.  At  this  scale  (0.1  Rjr),  the  models  are 
difficult  to  resolve,  thus  figure  5.3c  was  done  to  show  the  best  and  worst  agreement  at  this 
longitude. 

Figure  5.3c  is  at  an  altitude  of  15,000  km.  The  plot  on  the  left  is  done  at  the  same  sole  as 
figures  5.3a  and  5.3b.  The  two  boxes  drawn  on  this  plot  show  the  scale  of  the  two  plots  at 
the  right.  The  upper  panel  has  a  scale  10  times  larger  than  the  lower  panel.  L  values  were 
only  obtained  at  5°  increments,  therefore  the  close-up  figures  only  show  four  points  rather 
than  lines  in  order  to  see  some  separation  between  the  points.  In  the  upper  plot, 
Barrndough  is  roughly  .035  RE  higher  than  IGRF  and  Cain  is  roughly  .0025  RE  higher  than 
1GRF.  This  first  difference  is  on  the  lower  edge  of  detectability  by  CRRES  instruments. 
I  lowever,  if  data  is  binned  by  1/20  RE  as  has  been  discussed  (at  least  at  a  preliminary  survey 
level)  then  even  this  difference  is  negligible.  ‘Hie  second  difference  is  too  small  to  matter  to 
CRRES  at  any  level.  Note,  this  upper  plot  is  on  a  scale  an  order  of  magnitude  larger  than 
the  lower  plot,  thus  these  differences  arc  also  of  no  consequence  to  CRRES.  In  fact,  IGRF 
.md  Cain  plot  on  lop  of  each  other  in  this  lower  plot. 


5.3  L  DETERMINED  BY  IGRF  1985  FOR  1990 

Contour  plots  of  L  over  the  entire  latitude  and  longitude  range  are  shown  in  Figures  5.4a*c. 
Figure  5.4a  shows  the  temporal  changes  in  L-shell  as  found  from  IGRF  1985  evaluated  at 
350  km  for  1990  and  19S0.  There  is  no  significant  change  seen  here.  In  fact,  if  one  goes 
back  to  plots  generated  in  1970  by  Stassinopoulos,  one  sees  no  change  over  20  years. 
Figures  5.4b  and  5.4c  show  L-shell  contours  from  IGRF  1985  evaluated  for  1990  for  the 
other  four  primary  altitudes  considered  in  this  report  (850  km,  1500  km,  10,000  km,  and 
15,000  km).  These  are  included  as  a  general  reference  for  the  reader.  Finally,  Table  5.1 
shows  the  maximum  differences  between  the  models  for  1980  and  1990  at  six  altitudes  (0.0 
km  to  10,000  km  at  2000  km  increments).  These  dearly  show  that  Barraclough  is  furthest 
from  IGRF  and  Cain  is  closest  for  1990.  For  1980,  the  differences  are  lower  overall  (as  is 
anticipated),  however,  MAGSAT  and  Cain  agree  about  eqaliy  well  with  IGRF.  Note,  this  is 
not  surprising  since  Cain’s  model  was  based  on  this  particular  MAGSAT  model  then 
supplemented  it  with  ground-based  data  (see  Giapter  2). 
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LONCITUDE  (DECREES) 


Figure  5.4a.  Contours  of  L-shcll  as  determined  by  IGRF  1985  over  the  full  range  of  latitudes  and 
longitudes  at  350  km  for  1990  (top)  and  19f0  (bottom).  There  is  no  visible  difference  between  these 
epochs. 
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Figure  5.4b.  Contours  of  L-shcil  as  determined  by  IGRF 1985  for  1990  over  the  full  range  of  latitudes  and 
longitudes  at  850  km  (top)  and  1500  km  (bottom). 
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Figure  5.4c.  Contours  of  L-shcll  as  determined  by  IGRF  19S5  for  1990  over  the  full  range  of  latitudes  and 
longitudes  at  10,000  km  (top)  and  15,000  km  (bottom). 
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Table  5.1.  MAXIMUM  L-SHELL  DIFFERENCES  BETWEEN  ±20°  LATITUDE 
All  models  Compared  to  IGRF  1985  (ModeMGRF) 


Year 

Alt.  (xlO3  km) 

Barr  (xlO3) 

Magsat  (xlO-3) 

Gun  (xlO ') 

IV.SO 

0.0 

4.10 

1.32 

2.86 

10 

4.13 

1.08 

1.98 

4.0 

4.42 

1.13 

1.61 

6.0 

5.13 

1.03 

•0.57 

8.0 

4.98 

1.28 

1.04 

10.0 

4.73 

1.06 

-0.54 

1W0 

0.0 

19.89 

11.19 

3.26 

10 

20.31 

10.57 

2.10 

4.0 

2150 

10.40 

2.45 

6.0 

24.59 

11.27 

2.44 

ao 

25.99 

11.59 

16-1 

10.0 

2a77 

1130 

3.51 
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6,  COMPARISONS  OF  THE  SOFTWARE 


Another  aspect  of  the  analysts  of  these  models  is  that  of  the  software  which  uses  them.  A 
package  already  exists  for  CRRES  which  incorporates  the  IGRF  1985  model.  This  package 
(OF! ’RACE)  was  modified  to  include  the  other  three  models:  Barradough  1S75,  MAGSAT 
19S0,  and  Cain  1990.  The  MAGSAT  model  contains  50  coefficients  to  allow  for  more 
detailed  modelling  of  the  field.  However,  evaluation  of  this  many  terms  requires  an 
exorbitant  amount  of  computer  time.  Thus,  a  truncated  version  of  this  model  was  used  in 
this  analysis.  OPTRACE  was  further  modified  to  exclude  external  field  contributions.  'Hus 
routine  requires  more  CPU  time  than  another  code  also  in  use  at  GL,  HMIN.  This  second 
code  was  modified  to  give  the  same  results  as  those  initially  obtained  from  OPTRACE  (see 
Chapter  4).  It  was  this  code  which  was  used  to  do  the  B  field  and  L-shell  comparisons  in 
Chapters  3  and  5,  respectively. 

So  ft  wire  was  received  from  Dr,  J.  C  Cain  to  calculate  L-sbell  using  both  integration  of  the 
spherical  harmonic  coefficients  and  using  interpolation  tables.  Both  of  these  codes  from 
Ciin  were  originally  developed  in  the  early  1970s  by  Kluge  (see  Owpter  4).  In  order  to 
avoid  confusion,  the  integration  code  will  be  referred  to  as  Kluge/Cain  and  the  interpolation 
code  will  be  referred  to  as  Interpolation,  ’live  advantage  of  using  the  interpolation  tables  is 
one  of  greatly  reduced  CPU  time  in  the  processing.  However,  the  question  arises  as  to 
whether  or  not  the  increased  speed  sacrifices  accuracy.  In  comparisons  between  interpolated 
L  values  and  those  obtained  by  integration,  it  is  found  that  the  interpolated  values  are 
sufficiently  accurate  for  most  applications.  However,  a  problem  was  found  with  these  codes. 
Oiin  coded  a  value  for  the  dipole  moment  into  the  integration  code  rather  than  calculating  it 
from  the  coefficients  (see  Chapter  4).  The  value  he  selected  dates  back  to  the  mid-1960s  and 
is  now  out  of  date.  The  dipole  moment  and  the  position  of  the  offcet  dipole  are  changing 
rapidly  enough  that  a  noticeable  difference  is  seen  over  20  years  (see  Section  4.3.3.1).  'Ihu«, 
the  interpolated  values  are  not  as  accurate  as  they  could  be  since  the  tables  were  based  on 
the  results  from  Kluge/Cain.  Using  results  from  a  model  which  incorporates  the  current 
dipole  moment,  a  more  accurate  set  of  tables  could  be  obtained. 

A  similar  approach  to  that  used  with  the  magnetic  field  and  L-shell  comparisons  was  taken. 
To  begin  with,  contour  plots  of  the  differences  between  tire  codes  were  generated.  Then  line 
plots  to  high  latitude  at  selected  longitudes  were  made.  More  detailed  lincplots  within  the 
CRRES  latitude  range  were  then  considered.  After  reviewing  the  discrepancies  between  the 
codes,  their  efficiency  was  also  evaluated. 


6.1  OVERALL  COMPARISONS  OF  TI  IE  CODES 

Contour  plots  over  the  entire  latitude  and  longitude  ranges  of  the  differences  between 
Interpolated-OPTRACE/Cain,  Interpolated-OPTRACE/IGRF,  InterpoIated-KJuge/Cain,  and 
Klugc/Cain-OPTRACE/Cain  are  shown  in  Figures  6.1  and  6.2.  Figure  6.1a  shows  the  L-shell 
differences  for  these  four  cases  for  1990  and  350  km  altitude.  Three  of  the  plots  are  to  be 
used  in  evaluating  the  accuracy  of  the  Interpolation  code.  The  fourth  plot  (lower  right 
panel)  shows  the  agreement  between  Kluge/Cain  and  OPTRACE  using  the  same  model. 
Kluge/Cain  uses  a  transformation  to  inverse  coordinates  to  simplify  the  integration  (i.e.,  the 
integration  is  performed  over  a  straight  line  rather  than  a  curve  as  is  done  in  OPTRACE). 
However,  one  expects  the  agreement  here  to  be  quite  dose.  The  agreement  is  seen  to  be 
very  gocJ  in  all  four  plots.  Recall  that  in  the  model  comparisons,  the  differences  between 
Cain  ana  IGRF  were  very  small  and  for  the  most  part  negligibly  small.  This  is  apparent 
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-Shell  Difference  Contours  :  Year  -  1990.0  Attitude  —  350.00 


(s33H93a)  3anmvi 


(s33yo3a)  3aruiJL\n 


LONGITUDE  (DEGREES)  LONGITUDE  (DEGREES) 

Figure  6.1a.  Contour  plots  of  L-shcll  differences  in  Earth  radii  (R^  covering  the  full  latitude  and  longitude 
ranges  as  determined  form  the  various  software  packages:  Interpolated,  Kluge/Cain,  OPTRACE/Cain,  3r.d 
OPTRACE/IGRF.  The  comparisons  arc  done  for  an  altitude  of  350  km  (perigee)  at  1990.0. 
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Figure  6.1b.  Contour  plots  of  model  L-slicll  differences  in  Earth  radii  (R,^  covering  the  full  latitude  and 
longitude  ranges  as  determined  form  the  various  software  packages:  Interpolated,  Klugc/Cain, 
OPTRACE/Cain,  and  OPTRACE/1GKF.  The  comparisons  arc  done  for  an  altitude  of  10,000  km  at  1990.0. 
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here  as  (he  (op  two  panel*  which  show  the  difference*  between  Interpolated  and 
t 'ptrace/Cain  and  OPTRACE/1GRF,  respectively,  show  very  similar  difference  contours, 
the  differences  between  Interpolated  and  Khige/Grin  are  very  small  (most  of  the  contours 
are  0.  level  contours).  This  is  expected  since  the  Interpolation  tables  were  derived  from  the 
results  of  Kluge/Caln.  ‘thus,  the  differences  between  Kluge/Cain  and  OPTRACE/Ctnn  (lower 
right  panel)  are  nearly  the  same  as  those  seen  for  Interpolated  and  OPTRACE/Cain  (upper 
left  panel). 

These  same  properties  are  seen  in  Figure  6,1b  as  welt  Here,  comparisons  are  made  at  an 
altitude  of  10,000  km.  The  differences  arc  a  bit  smaller  than  at  350  km.  Again,  the  top  two 
plats  and  the  lower  right  plot  look  very  similar  due  to  the  close  agreement  between  the  Gain 
anJ  IGRF  models  ami  (he  close  agreement  between  Interpolated  and  Kluge/Cain.  Again, 
the  smallest  differences  are  found  between  Interpolated  and  Kluge/Cain  (Tower  left  plot). 
However,  when  the  discrepancy  between  the  dipole  moments  arose,  these  results  changed  a 
httlc  bit.  Re-evaluaiing  Kluge/Cain  with  the  updated  moment  led  to  better  agreement 
between  KIuge/Qin  and  OPTRACE/Cain  than  is  seen  for  Interpolated  and  Khige/Giitt. 
Note,  one  is  not  able  to  similarly  update  the  interpolation  tables,  they  must  be  regenerated 
from  results  which  were  obtained  with  the  current  moment,  lit  us,  in  Figures  6.2a  and  6.2b 
(350  km  and  10,0000  km,  respectively),  the  picture  has  changed  somewhat  from  Figures  6.1a 
,md  6.1b.  One  now  sees  the  best  agreement  between  Kluge/Cain  and  OPTRACE/Cain 
(lower  right  panel)  and  the  other  three  plots  have  very  similar  difference  contour 
distributions.  Note,  the  differences  have  not  increased. 


f2  COMPARISONS  AT  SPECIFIED  LONGITUDES 

liter  high  latitude  line  plots  taken  at  lour  longitudes  over  the  globe  (20®,  150®,  200®,  and  290° 
east  longitude)  show  the  results  from  Interpolated,  Kluge/Cain,  OPTRACE/Cain, 
OPTRACE/IGRF  (Figure  6.3).  At  this  resolution,  no  difference  between  the  four  routines 
can  be  seen. 

A  plot  done  over  the  ±20°  latitude  range  shows  the  offset  of  the  Interpolated  results  from 
the  those  of  the  other  three  (note  the  current  moment  was  used  in  Kluge/Gain  here)  in 
Figure  6.4.  The  Interpolated  results  are  represented  by  the  solid  litre  at  higher  L  values  than 
the  other  three  (which  appear  to  be  one  line,  the  dashed,  dotted,  and  dot-dash  lines  all  plot 
over  each  other).  When  this  was  first  done  with  the  old  moment,  Interpolated  and 
Kluge/Cain  plotted  over  each  other  and  were  offset  from  the  OPTRACE  results  to  higher 
values  of  L  This  is  what  is  expected  since  a  higher  dipole  moment  leads  to  higher  values  of 
L 

To  separate  the  models  a  little  better,  a  plot  with  details  on  a  much  smaller  scale  was  done, 
Figure  6.5.  Here  a  line  plot  like  that  shown  in  figure  6.4  was  done  at  80®  east  longitude  und 
15,000  km  altitude  for  1990.  TWo  boxes  are  drawn  on  the  plot  on  the  left  to  show  the 
regions  which  are  detailed  on  the  right.  Figure  6.5  was  produced  using  tiie  old  dipole 
moment.  Thus,  in  both  details,  Interpolated  and  Kluge/Crin  are  very  near  each  other  and 
tire  shifted  to  a  higher  L  than  the  OPTRACE  results  which  are  also  near  each  other.  Note, 
the  upper  panel  has  a  scale  which  is  twice  as  large  as  that  of  the  lower  plot.  The  differences 
are  small,  but  not  negligible.  In  the  upper  plot,  Interpolated  and  Kluge/Cuin  are  about  .04 
Ri;  different.  In  the  lower  plot,  the  difference  is  roughly  .025  for  Interpolated  and  .03  for 
Kluge/Cbin. 
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Figure  6.4.  Line  plots  similar  to  those  in  Figure  6.3,  but  over  a  smaller  latitude  range.  Even  at  this 
resolution,  the  differences  between  the  models  are  very  small.  Here,  the  current  value  for  the  dipole 
moment  has  been  used,  thus,  only  Interpolated  results  are  offset  from  the  other  three  routines.  The  top 
panel  shows  the  differences  at  350  km  and  SO*  cast  longitude  for  1990.0.  The  bottom  panel  is  the  same 
plot  done  fer  epoch  10,000  km.  The  offset  is  larger  for  the  higher  altitude,  but  KJuge/Cain, 
OPTRACE/Cain,  and  OPTRACE/IGRF  still  agree  very  well. 
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were  found  using  lire  old  dipole  moment.  Note,  how  well  Interpolated  and 
Kluge/Gun  agree  and  how  far  they  arc  offset  from  OPTRACE/Giin  and 
OPTRACE/IGRF. 


Using  (he  current  dipole  moment  in  Kltigc/Cain  brings  it  into  closer  agreement  with  the 
OPTRACE  results  as  seen  in  Figure  6.6.  In  the  left  panel  of  figure  6.6,  Kluge/Cain  now 
coincides  with  the  OPTRACE  results.  The  upper  detail  shows  a  shift  of  about  .03  closer 
to  the  OPTRACE  points.  The  lower  detail  has  a  slightly  larger  shift  than  that,  placing 
Kluge/Cain  a  little  bit  lower  than  the  OPTRACE  points.  Clearly,  the  different  dipole 
moment  makes  a  significant  difference  in  the  L*shell  results  (see  Table  6.1) 


6.3  TIMING  COMPARISONS 


6.3.1  Optimizing  the  Code 

At  the  start  of  this  project,  OP'IRACE  was  modified  to  perform  internal  field  calculations 
only  without  any  externa]  contributions.  It  was  further  modified  to  include  the  other  three 
models  (Barraclough,  MAGSAT,  and  Cain).  The  MAGSAT  model  was  comprised  of  n-50 
terms,  which  were  very  time  consuming  to  evaluate.  Thus,  MAGSAT  was  used  in  three 
different  subroutines  for  the  modified  OPTRACE  which  allowed  the  user  to  select 
interactively  whether  to  use  the  full  blown  n-50  set  or  one  of  the  two  smaller  sets,  one 
truncated  at  n-15  and  the  other  at  n-20.  Throughout  all  of  the  modifications,  tests  were 
run  to  check  for  accuracy  and  efficiency.  Up  until  the  full  field  comparisons  were  run  for 
each  of  the  models  (using  the  n-15  MAGSAT  version),  short  tests  were  run  on  a  selected 
few  points  which  spanned  the  latitude,  longitude  and  altitude  range.  It  was  found  that  the 
differences  were  small  between  MAGSAT  n-15  and  n-50.  Even  though  the  detail  of  the 
higher  order  model  was  lost,  the  time  saved  justified  it  (Table  6.2). 

Another  factor  of  interest  when  discussing  efficiency  is  that  of  the  number  of  segments  (or 
steps)  used  in  integrating  over  the  field  line.  The  more  segments  used,  the  better  the 
accuracy  will  be.  However,  this  can  also  become  time  consuming,  so  one  needs  to  sacrifice 
some  of  the  accuracy  for  the  sake  of  efficiency  without  compromising  the  integrity  of  the 
results.  For  the  purpose  of  this  report,  50  segments  were  found  to  be  sufficient  (Table  6.3). 


6.3.2  Comparisons  with  Cain’s  Routines 

Having  resolved  these  questions,  the  full  field  comparisons  were  done  using  OPTRACE  with 
50  segments  (Table  4).  Clearly,  these  were  time  consuming  computations.  In  comparing  the 
efficiency  of  OPTRACE  with  the  two  codes  sent  by  Cain,  it  became  necessary  to  remove  all 
of  the  extraneous  calculations  done  in  OPTRACE  so  that  only  L  was  being  determined. 
This  would  have  proved  to  be  an  intricate  task,  so  to  save  a  good  bit  of  effort,  HMIN  was 
modified  for  the  purpose  of  simply  finding  L  HMIN  uses  essentially  the  same  routines  as 
OP'IRACE  and  by  replacing  two  subroutines,  it  was  made  into  the  same  algorithm  as  that 
in  OPTRACE  for  consistency  sake.  There  is  a  difference  regarding  the  length  along  the 
field  line  each  of  these  routines  traces.  OPTRACE  traces  to  tire  100  km  intercept  whereas 
HMIN  goes  only  to  the  conjugate  point  (see  Chapter  4).  Thus,  HMIN  is  the  faster  *  .f  the 
two  routines. 

For  the  full  fields  obtained  for  1965,  the  modified  Hmin  routine  was  used  to  obtain  L  The 
CP  time  required  to  do  so  was  then  compared  to  the  Kluge/Cain  integration  method  and  to 
th  Interpolation  tables  (Table  6.5).  For  internal  fields  only,  the  interpolation  is  the  most 
efficient.  However,  to  really  take  advantage  of  this,  one  needs  to  update  the  tables  with  the 
current  dipole  moment.  The  integration  sent  by  Cain  is  also  significantly  faster  than  the 
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Tabic  (\1.  Changes  in  the  Dipole  Moment  and  Position  as  Gilculatcd  by  Barrnclnugli  1975, 
MACtSAT  1980,  IGRF  1985.  and  Ciin  1990  Extrapolated  to  1990.5. 1975.5,  and  1965.5. 


Dipole  AUributci 

Year 

Moment 

1990.5 

1975.5 

1965.5 

Offset  from  Center 

1990.5 

(in  UB) 

1975.5 

1965.5 

latitude 

1990.5 

1975.5 

1965.5 

Longitude 

1990.5 

1975.5 

1965.5 

Darraclough  Sfagiat 

0.30155  0.30270 

0.30687  0.30699 

0.30975  0.30987 


0.0S0Q6  0.081  14 

0.07460  0.074S6 

0.07118  0.07075 

13.6  21.8 

19.9  19.8 

17.2  18.3 

147,2  146.3 

147.9  147.7 

148.5  148.7 


IGRFJ9S5 

Cam 

0.30285 

0.30299 

0.30703 

0.30706 

0.309S3 

0.309SQ 

0.08100 

0.08135 

0.0750-1 

0.07415 

0.07114 

0.06962 

21.3 

21.2 

20.0 

20.1 

19.1 

19.3 

146.1 

146.0 

147.7 

148.0 

148.9 

149.4 
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umie  0.2.  Inferences  in  L  values  due  If?  using  more  coefficients  in  the  Macsat  model. 
Mure  coefficients  lead  to  more  detail,  particularly  at  low  altitudes.  However,  based  on  the 
comparisons  below,  the  improvement  did  not  justify  the  additional  computation  time 
•Ci|»nCu  uy  die  greater  number  of  coefficients  (9.40  times  longer). 


GUT 

GLON 

A  IT 

L-SHEL(lS) 

L-SIIEL(SO) 

U(IS-SO) 

»«, 

36. 

a 

28.13959 

28.14000 

-.00041 

Til! 

72. 

a 

5.40458 

5.40438 

JWOM 

JO. 

m 

a 

164736 

1J54930 

.00006 

30. 

144. 

a 

1.10648 

1.10647 

.00001 

to. 

mu. 

0. 

.98319 

.98320 

-.00001 

0. 

216. 

a 

.98005 

.97992 

.00013 

•20, 

252. 

0. 

13)5520 

1.05523 

•00003 

-to. 

286 

a 

1.27467 

1.27473 

-.0X06 

•00. 

324. 

a 

2.21630 

121627 

.00003 

■SO, 

360. 

a 

656247 

6.70409 

-.14162 

85. 

36. 

500. 

2952115 

2952111 

.00004 

70, 

72. 

500, 

5.64643 

5.64677 

.00005 

JO. 

106 

500. 

1.96966 

1.96965 

.00001 

30. 

144. 

500. 

1.19406 

1.19406 

... 

10. 

laa 

500. 

1.05940 

1.05940 

... 

a 

216. 

500. 

1.05868 

1.05867 

.00001 

•20. 

252. 

500. 

1.13369 

1.13370 

-.00001 

•40. 

288. 

500. 

157512 

157512 

... 

00. 

324. 

500. 

2.41738 

141739 

>.00001 

■m 

360. 

500. 

7.43303 

7.43300 

3)0003 

ss. 

36. 

1000. 

3053374 

3053370 

.00004 

7a 

72. 

1000. 

5.90096 

5.90094 

.00002 

so. 

108 

1000. 

2.09167 

109166 

.00001 

30. 

144. 

1000. 

1.28178 

1.28178 

ta 

180, 

1000. 

1.13587 

1.13587 

••• 

0. 

216. 

1000. 

1.13705 

1.13705 

••• 

■20. 

251 

1000. 

1.21287 

1.21287 

-to. 

288, 

1000. 

1.47585 

1.47585 

... 

•60. 

324. 

1000. 

161466 

16)466 

••• 

•SO. 

360. 

1000. 

851703 

851702 

.00001 

85. 

36. 

5000. 

41.04362 

41.04362 

••• 

70. 

71 

5000. 

8.19541 

8.19541 

... 

50. 

108. 

5000. 

3.09190 

3.09190 

... 

30. 

144. 

3000. 

1.98775 

1.98775 

... 

10. 

180. 

5000. 

1.75498 

1.75498 

... 

0. 

216. 

5000. 

1.76321 

1.76321 

•20. 

251 

5000. 

1.85660 

1.85660 

... 

•40. 

288. 

5000. 

128891 

128891 

... 

•60. 

324. 

5000. 

4.16962 

4.16962 

•80. 

360. 

5000. 

15.46061 

15.46061 

... 

85. 

36. 

10000. 

54.94462 

54.94462 

70. 

71 

10000. 

11.28203 

11.28205 

... 

50. 

108. 

10000. 

456666 

456666 

... 

30. 

144. 

10000. 

187586 

187586 

... 

10. 

180. 

10000. 

253910 

153910 

... 

0. 

216. 

10000. 

154807 

154807 

... 

•20. 

252. 

10000. 

166963 

166963 

... 

•40. 

288. 

10000. 

350894 

350894 

M. 

•60. 

324. 

10005. 

6.09631 

6.09631 

... 

•80. 

360. 

10000. 

24.41243 

24.41243 

... 

85. 

36. 

15000. 

69.10153 

69.10153 

... 

70. 

71 

15000. 

14.43114 

14.43114 

50. 

108. 

15000. 

5.65149 

5.65149 

... 

30. 

144. 

15000. 

3.76838 

3.76838 

... 

10. 

180. 

15000. 

352745 

352745 

... 

0. 

216. 

15000. 

353417 

353417 

... 

•20. 

251 

15000. 

3.48378 

3.48378 

... 

•40, 

288. 

15000. 

452936 

452936 

... 

•60. 

324. 

15000. 

8.02087 

8.02087 

... 

-80. 

360. 

15000. 

3355608 

3355608 

— 
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fable  6.3.  Differences  in  L  values  due  to  using  more  segments  in  the  integration  of  the  field 
line.  More  segments  lends  to  greater  accuracy,  homer  much  more  CPU  time  is  then 
required  to  do  the  evaluation  (1.77  times  ionger).  The  improvement  of  100  segments  over  50 
segments  was  not  deemed  sufficient  to  warrant  the  additional  time  required  for  the  full  held 
•.'unpansons  for  this  report. 

CUT  GLON  ALT  IGRF(SEGS=SO)  IGRF(SEGS-IOO)  DIFFERENCE 


85 

36. 

4 

28.02246 

28.05233 

-.02965 

7U. 

72. 

4 

5.39173 

5495418 

s004IS 

M 

108. 

0. 

144718 

144766 

•.00O»8 

3d 

144. 

4 

1.10581 

1.10583 

-JM0Q2 

to 

ISO. 

4 

56299 

.96300 

-.00001 

0. 

216. 

4 

.96014 

.96014 

.00000 

•a 

251 

0. 

105618 

1.05619 

-.00001 

2*8. 

0. 

1.27461 

1.27470 

-.00009 

•<4 

324. 

4 

121627 

121781 

-.00154 

*84 

360. 

4 

645010 

656058 

-.01018 

85. 

36. 

500. 

29.19632 

29.22007 

-.02175 

74 

72. 

500. 

5.63183 

5.63594 

-.00111 

50. 

108. 

500. 

1.96743 

1.96838 

-.00095 

30. 

144. 

500. 

1.19340 

1.19346 

-.00006 

14 

ISO. 

500. 

1.05930 

1.05931 

-.00001 

0. 

216. 

504 

1.05877 

145878 

-.00001 

■2d 

252. 

500. 

1.13447 

1.13450 

-.00003 

-to 

288. 

500. 

147503 

147518 

-.00015 

to 

324. 

500. 

141701 

141791 

-.00090 

l >0 

360. 

500. 

7.42500 

7.43425 

-.00925 

85. 

36. 

1000. 

30.41567 

30.42866 

-.01299 

7il 

71 

1000. 

548853 

5.89065 

-.00232 

54 

108. 

1004 

108936 

109028 

-4)0092 

30 

144. 

1000. 

1.28112 

1.28121 

-.00009 

10. 

160. 

1000. 

1.13582 

1.13584 

-.00002 

u. 

216. 

1004 

1.13716 

1.13718 

-.00002 

•20. 

251 

1000, 

1.21351 

1.21355 

-.00004 

-14 

288. 

1000. 

1.47574 

1.47620 

-.00046 

*60. 

324. 

1000. 

161405 

161580 

-.00175 

-SO. 

360. 

1000. 

840901 

841489 

-.00588 

85. 

36. 

5000. 

4487692 

4048696 

-.01006 

70. 

71 

5000. 

8.18040 

8.18469 

-.00429 

50. 

108. 

5000. 

3.06894 

3.06966 

-.00092 

30. 

144. 

5000. 

1.96704 

1.96756 

-.00052 

10 

180. 

5000. 

1.75502 

1.75535 

-.00033 

0. 

216. 

5000. 

1.76341 

1.76358 

-.00017 

*20. 

251 

5000. 

145671 

145728 

-.00057 

-10. 

288. 

5000. 

128654 

128959 

-.00105 

-60. 

324. 

5000. 

4.16793 

4.16926 

-.00133 

•SO. 

360. 

5000. 

15.44047 

15.44737 

-.00690 

85. 

36. 

10000. 

54.72477 

54.73180 

-.00703 

70. 

72. 

10000. 

11.26230 

11.26473 

-.00243 

50. 

108. 

10000. 

446294 

446490 

-.00196 

30. 

144. 

10000. 

187506 

187641 

-.00135 

10. 

180. 

10000. 

153915 

153952 

-.00037 

0. 

216. 

10000. 

154827 

154857 

-.00030 

•20. 

251 

10000. 

166947 

167000 

-.00053 

-10. 

288. 

10000. 

340815 

340901 

-.00086 

-60. 

324. 

10000. 

6.09358 

6.09528 

-.00170 

-80. 

360. 

10000. 

2448652 

2449253 

-.00601 

85. 

36. 

15000. 

6842715 

68.83137 

-.00422 

70. 

71 

15000. 

14.40627 

14.40871 

-.00244 

54 

108. 

15000. 

5.64692 

5.64828 

-.00136 

30. 

144. 

15000. 

3.76747 

3.76839 

-.00092 

10. 

180. 

15000. 

342752 

342827 

-.00075 

0. 

216. 

15000. 

343439 

343473 

-.00034 

•20. 

251 

15000. 

3.48343 

3.48406 

-.00065 

-14 

288. 

15000. 

442814 

442962 

-.00168 

•60. 

324. 

15000. 

8.01706 

8.01896 

-.00190 

-80. 

360. 

15000. 

3341976 

3342451 

-.00475 
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'I'i.ble  6.4.  Ratio*  of  the  required  CPU  time  for  die  foil  field  L-shell  evaluations  using  the 
four  models  in  OPTRACE.  L  determined  for  every  5°  in  latitude  over  ±75°  and  for  every 
10°  in  longitude.  Various  altitudes  were  grouped  together  for  the  computations  as  indicated. 
The  average  computation  time  for  altitudes  between  0  km  and  10,000  km  using  IGRF  19.85 
extrapolated  to  196$  was  used  as  the  unit  time. 


YEAR 

Barracleugk 

Magiat 

IGRF  1985 

Cain 

50Gkm-1000km 

I960 

1.279 

1.874 

0.951 

1.865 

(step-500km) 

19S4 

1.290 

1.871 

0.960 

1.878 

19S8 

1.289 

1.857 

0.949 

1.854 

Qkm-10, 000km 

19S0 

1.380 

1011 

2.014 

(step"* 2000km) 

1982 

1.379 

2.001 

1.020 

2.013 

19S4 

1.3S3 

2.002 

1.024 

1010 

1986 

1.376 

2.016 

1.024 

1014 

1988 

1.378 

1009 

1.017 

1003 

1990 

1.373 

1014 

1.018 

100-1 

I5,000km-40,000km 

1980 

2.204 

3.220 

1.641 

3.219 

(step-5000km) 

1984 

2.211 

3.229 

1.648 

3.239 

1988 

2.210 

3.227 

1.6-15 

3.227 

Table  6.5.  Ratios  of  the  CPU  lime  required  by  HMIN,  Kluge/Cain  Integration,  and 
Interpolation.  These  arc  full  field  evaluations  of  L  taken  every  5°  over  ±75°  in  latitude  and 
over  ever>*  10°  in  longitude.  ’Hie  computations  were  done  for  various  altitudes,  grouped  as 
indicated.  All  evaluations  were  for  1985.  The  average  computation  time  for  altitudes 
between  0km  and  10,000km  using  IGRF  1985  extrapolated  to  1988  was  used  as  the  unit  time. 


ALTS 

500km— 1000km 
(»tep=500) 

2000km-8000km 

($tep=2000) 

I0,000km-J0,000km 

(stcp=5000) 

Interpolated 

0.003 

0.003 

0.003 

Integrated 

0.052 

0.042 

0.044 

HMIN(Cain-50) 

0.812 

0.817 

1.443 

HMIN(IGRF-50) 

0.411 

0.411 

0.733 

I-IMIN(Cain-lOO) 

1.305 

1.268 

2.297 

HMIN(IGRF-IOO) 

0.663 

0.637 

1.138 
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technique  used  in  OPTRACE  or  HMIN.  This  is  due  to  the  use  cJ  inverse  coordinates 
centered  on  the  dipole  so  that  one  integrates  over  a  straight  line  rather  than  a  curve  ns  in 
OPTRACE, 
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